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ABSTRACT 


Deep  reflections  from  horizons  within  the  Earth 
have  been  recorded,  at  near  vertical  incidence,  on 
seismograms  along  three  lines  in  southwestern  Manitoba. 

The  data  were  digitally  recorded  on  magnetic  tape  recorders 
with  a  dynamic  range  of  84  decibels.  Power  spectral  esti¬ 
mates  indicate  that  the  reflected  energy  is  contained  in 
the  frequency  band  of  10  to  30  hertz.  Through  the 
application  of  frequency  filtering  and  common  depth  point 
stacking,  the  primary  reflections  were  enhanced  relative 
to  the  noise. 

The  methods  of  N-th  root  stacking,  velocity  filter¬ 
ing  and  maximum  likelihood  filtering  are  discussed.  These 
techniques  are  applied  to  the  reflection  data  and  the 
results  are  compared.  Velocity  filtering  and  N-th  root 
stacking  were  both  found  to  be  effective  in  enhancing 
coherent  signal,  with  very  little  success  being  obtained 
from  a  three  channel  maximum  likelihood  filter.  This  was 
the  first  time  that  the  N-th  root  stack  had  been  used  for 
enhancing  reflection  records.  Main  reflected  energy  was 
located  at  depths  of  about  19  (6.6  seconds),  24  (8.1 
seconds)  and  37  (12.5  seconds)  kilometers.  Also,  the  N-th 
root  stack  was  found  to  be  an  excellent  filter  for  deter¬ 
mining  the  phase  velocity  of  the  reflected  events. 
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CHAPTER  1 
INTRODUCTION 


1.1.  Understanding  the  Crust 

The  crust  is  comprised  of  the  uppermost  layers  of 
the  Earth.  This  outer  rind  constitutes  about  one  percent 
of  the  Earth's  volume,  and  only  0.5  percent  of  the  Earth's 
mass  (Jacobs  et  al .,  1974).  It  has  a  mean  thickness  of 
about  10  kilometers  under  the  oceans  and  40  kilometers 
under  the  continents,  and  consists  mainly  of  oxygen,  sili¬ 
con,  aluminum  and  iron.  Separated  from  the  crust  by  the 
Mohorovicic  discontinuity,  lies  the  mantle,  and  beneath 
this  is  the  core  (figure  1.1).  However,  to  fully  under¬ 
stand  these  deeper  regions,  in  particular  the  mantle,  one 
must  completely  understand  the  crust. 

Explosion  seismology,  and  in  particular  the  seismic 
reflection  technique,  has  proven  to  be  the  most  valuable 
method  for  investigating  deep  crustal  structure.  Both 
reflection  and  refraction  methods  have  produced  much  infor¬ 
mation  about  the  lower  crust  and  upper  mantle,  but  the 
reflection  technique  has  many  inherent  advantages.  First, 
the  wavelengths  of  the  reflected  waves  are  many  times 
shorter  than  those  of  refracted  waves  yielding  a  larger 
resolving  power  of  the  method,  and  thus  increasing  the  pre¬ 
cision  in  detailing  fine  structure  (Kanasewich  and  Cumming, 
1965).  Spectral  analysis  of  reflected  waves  provides 
information  about  the  reflective  process.  Also,  the 
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Figure  1.1:  Depths  to  the  major  discontinuities  of  the 
Earth 


The  dotted  region  indicates  complex  structure  in  the  upper  mantle. 
A  radius  of  6371  kilometers  is  the  result  of  assuming  the  Earth  to 
be  a  perfect  sphere  with  its  true  mass  and  volume. 


method  allows  the  determination  of  the  average  velocity 
down  to  any  discontinuity,  and  the  detection  of  all  low 
velocity  zones. 

Through  the  use  of  digital  recorders  and  high  speed 
digital  computers,  seismic  reflection  data  may  be  processed 
through  various  signal  enhancement  techniques  very  quickly 
and  inexpensively.  Thus,  seismic  sections  exhibiting  the 
very  fine  detail  of  deep  structures  within  the  Earth  may 
be  produced.  These  sections  may  then  be  correlated  with 
other  geophysical  information  obtained  from  refraction, 
magnetic  and  gravity  surveys,  to  produce  large  scale  pic¬ 
tures  of  the  Earth's  crust. 

1.2.  History  of  the  Reflection  Technique 

The  development  of  the  seismic  reflection  method 
began  in  1917  at  the  University  of  Oklahoma,  in  the  hopes 
of  developing  a  new  technique  to  locate  oil  reserves 
(Schriever,  1952).  Mr.  J.C.  Karcher  and  Dr.  W.P.  Haseman 
met  at  the  university,  and  discussed  the  advantages  and 
disadvantages  of  using  reflected  seismic  waves  to  locate 
oil  bearing  structures.  This  idea  was  left  alone  until 
1919,  when  they  again  met,  leading  to  the  design  and  con¬ 
struction  of  equipment  for  the  recording  of  seismic 
reflections.  The  world's  first  exploration  record  was 
made  in  1919  by  Dr.  Karcher  in  a  rock  quarry  near  Washing¬ 
ton  D.C.,  using  dynamite  as  the  source.  This  was  followed 
by  the  first  field  tests  on  June  4,  1921.  Drs.  J.C. 
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Karcher,  W.P.  Haseman,  I.  Perinne,  and  Mr.  W.C.  Kite  com¬ 
prised  this  first  field  crew.  To  this  day,  the  reflection 
technique  is  still  the  most  useful  petroleum  prospecting 
technique  (Dobrin,  1976). 

Deep  crustal  reflections  were  first  recorded  by 
Junger  (1951),  on  a  seismic  section  shot  in  Big  Horn  County, 
Montana.  He  had  observed  that  a  charge  of  11.4  kilograms 
would  create  a  long  lasting,  low  frequency  disturbance. 

To  determine  the  duration  of  this  pulse,  the  recording 
camera  was  allowed  to  run  for  ten  seconds.  A  distinct 
reflection  was  observed  at  8.5  seconds  on  the  record.  This 
experiment  was  repeated  in  adjacent  areas,  with  reflections 
being  observed  at  7.0  to  7.5  seconds.  Junger  effectively 
argued  that  these  reflections  must  be  reflections  from  near 
horizontal  layers  beneath  the  basement,  and  not  multiple 
reflections  from  shallower  horizons.  Other  reports  concern¬ 
ing  the  early  use  of  the  reflection  technique  for  studying 
the  lower  crust  may  be  found  in  Steinhart  and  Meyer  (1961). 
James  and  Steinhart  (1966)  followed  this  up  with  a  review 
of  seismic  reflection  studies  of  the  lower  crust  from  all 
parts  of  the  world,  for  the  period  of  1960  to  1965.  Some 
of  these  will  be  mentioned  in  the  following  pages. 

Since  the  work  of  Mintrop  (1949),  there  have  been 
many  reports  of  deep  crustal  reflections  in  Europe.  In 
these  early  reports,  the  criterion  for  an  event  to  be  con¬ 
sidered  a  reflection,  was  that  it  had  to  have  a  high 
apparent  velocity.  Many  such  records  obtained  from  France 
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and  Hungary  during  the  1950's  exhibited  such  features. 
However,  many  of  these  European  records  contain  only  poor 
quality  reflections.  This  is  the  case  with  Bath  and 
Tryggvason  (1962),  who  report  on  the  first  seismic  investi¬ 
gations  of  deep  crustal  structure  in  Fennoscandia.  Their 
project  was  to  investigate  near  vertical  incidence  reflec¬ 
tions  from  crustal  discontinuities.  The  sources  of  energy 
were  quarry  blasts  in  iron  ore  mines,  and  the  information 
was  recorded  with  refraction  equipment.  They  obtained 
depths  to  the  Conrad  and  Mohorovicic  discontinuities  of  19 
and  34  kilometers,  but  their  reflections  were  erratic  and 
very  weak.  In  Greece,  Papazachos  et  al.  (1966),  reported 
on  the  crustal  structure  near  Athens.  Through  the  analysis 
of  reflections  obtained  from  near  earthquakes,  he  produced 
a  three  layered  model  for  the  Earth's  crust.  Again,  though, 
his  results  were  not  very  well  substantiated. 

The  application  of  the  reflection  technique  for  deep 
crustal  studies  became  popular  in  Germany  when  Dohr  and 
Schulz  recorded  their  first  profiles  in  1954.  In  the  eval¬ 
uation  of  the  data,  Dohr  applied  statistical  methods  to 
the  deep  reflections  (Dohr  and  Fuchs,  1967).  He  found  that 
on  many  occasions  late  reflections  were  observable,  but 
that  it  was  impossible  to  correlate  them  over  more  than  400 
meters.  Using  several  hundred  seismograms,  histograms  of 
the  number  of  deep  reflections,  in  time  intervals  of  0.1 
to  0.2  seconds,  for  particular  regions  of  about  100  to  300 
square  kilometers  were  plotted  against  arrival  time.  Some 
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of  these  histograms  exhibited  one  or  more  peaks  which 
Do hr  interpreted  as  reflected  events  from  intracrustal 
discontinuities  and  the  Mohorovicic  discontinuity. 

Liebscher  (James  and  Steinhart,  1966)  later  applied  these 
statistical  techniques  to  a  number  of  areas  in  southern 
Germany.  The  histograms  constructed  for  these  areas  also 
exhibited  many  peaks.  Liebscher  applied  a  chi-square  test 
to  the  peaks,  and  found  them  significant  at  the  99  percent 
level.  The  statistical  results  of  these  tests  suggest 
that  he  has  indeed  recorded  primary  reflections  from 
within  the  Earth. 

Do  hr  and  Fuchs  (1967)  also  discuss  the  statistical 
nature  of  deep  seismic  reflections  as  performed  on  the 
several  thousands  of  seismograms  collected  up  to 
1965.  By  performing  chi-square  tests  on  many  of  the  histo¬ 
grams,  the  computed  value  of  chi-square  for  the  peaks  was 
found  to  be  significant  at  the  95  percent  level.  This 
indicates  that  the  variations  in  the  distributions  are 
larger  than  what  would  result  from  random  variations.  In 
most  cases,  there  were  three  main  peaks  observed.  The  width 
of  the  peaks  (0.5  to  1.3  seconds)  was  suggested  to  be  the 
result  of  the  boundary  between  any  two  of  the  layers  being 
a  transition  zone,  which  Dohr  and  Fuchs  estimated  to  be 
1.5  to  4.0  kilometers  in  thickness.  They  also  noted  that, 
in  recording  deep  reflections,  better  results  were 
obtained  when  the  recording  arrays  consisted  of  groups  of 
geophones  with  higher  frequencies  than  the  corresponding 
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refraction  equipment. 

To  verify  the  conclusions  drawn  from  the  statistical 
method,  a  wide  angle  reflection  survey  was  conducted  in 
southern  Germany  in  1964,  utilizing  the  common  depth  point 
technique.  Meisnner  (1966)  reported  on  the  initial  results 
of  this  survey.  Four  different  reflecting  horizons,  the 
basement,  the  Conrad  discontinuity,  a  sub-Conrad  discon¬ 
tinuity,  and  the  Mohorovicic  (Moho)  discontinuity  were 
noted,  with  the  Moho  being  the  strongest  reflector.  The 
travel  time  curves  were  not  exactly  hyperbolic,  and  he 
explained  this  as  the  result  of  the  lack  of  sharp  boundaries 
between  the  horizons,  with  a  gradient  zone  2.0 ±.7  kms.  thick. 

In  1968,  seismic  reflection  measurements  were  taken 
along  a  17  kilometer  long  profile  in  order  to  investigate 
the  structure  of  the  Ries  Crater  in  southern  Germany.  The 
seismic  signals  were  recorded  digitally  for  up  to  14  seconds, 
with  the  digitizing  interval  being  4  milliseconds.  A  total 
spread  length  of  1380  meters  was  used,  and  each  subsurface 
point  had  300  percent  coverage.  After  stacking,  frequency 
and  optimum  filtering,  the  main  reflection  zones  were  noted 
at  3  to  4  seconds,  6  to  7  seconds,  and  9  to  10  seconds 
(Angenheister  and  Pohl,  1971).  This  last  reflection 
corresponds  to  the  Mohorovicic  discontinuity. 

Two  six  kilometer  profiles  were  conducted  in  the 
Rhinegraben  area  by  Dohr  as  part  of  the  Upper  Mantle  Pro¬ 
ject  in  1971  (Dohr  and  Meissner,  1975).  Then,  in  1973,  a 
deep  crustal  project  was  performed  across  a  large  fault 
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zone  separating  the  Rhenish  Massif  and  the  Saar-Nahe 
trough  in  southern  Germany.  A  21  kilometer  reflection 
profile  across  the  fault  zone  was  recorded.  Geophones 
were  organized  into  a  36-unit  arrangement,  and  a  300  per¬ 
cent  subsurface  coverage  was  obtained.  Charges  were  up  to 
100  kilograms.  The  spread  length  was  2.4  kilometers,  and 
24  traces  were  recorded  per  shot.  Many  high  frequency 
reflections  from  the  sediments  were  recorded,  along  with 
reflections  from  the  Conrad,  a  sub-Conrad,  and  the  Mohoro- 
vicic  discontinuities.  The  fault  zone  was  observed  down 
to  the  crust-mantle  transition  zone,  separating  two  large 
blocks  of  variable  dipping  crustal  layers  (Bartelsen  et  al. , 
1976). 

Deep  crustal  studies  using  reflected  waves  in  the 
Soviet  Union  were  first  proposed  by  G.A.  Gamburtsev  in 
1939  (Zverev,  1967).  From  1948  to  1950,  he  directed  the 
first  experiments  for  the  recording  of  reflected  waves  from 
horizons  within  the  lower  crust.  In  1956,  the  deep  seismic 
sounding  method  was  adopted  for  large  scale  seismic  pro¬ 
jects.  By  1961,  more  than  5000  kilometers  on  land  and 
11,000  kilometers  on  sea  of  seismic  reflection  profiles, 
had  been  recorded. 

From  the  period  of  1962  to  1965,  many  Soviet  geo¬ 
physicists  studied  the  crust  beneath  the  Black  and  Caspian 
Seas  (James  and  Steinhart,  1966).  They  noted  that  both  seas 
did  not  have  the  granitic  layer  (P-wave  velocity  of  6.0 
kilometers  per  second)  beneath  them,  whereas  the 


surrounding  regions  did.  The  crust  beneath  the  Black  Sea 
appeared  to  be  basaltic  and  at  most  20  kilometers  in 
thickness.  Beneath  the  Caspian  Sea,  the  crust  was  found 
to  be  more  acidic. 

Sollogub  (1973)  presents  the  results  of  deep  seismic 
sounding  in  central  and  southeastern  Europe.  The  study  is 
based  on  approximately  13,000  kilometers  of  reflection 
seismic  profiles.  Crustal  thicknesses  were  found  to  vary 
from  20  to  65  kilometers.  He  suggest.s  that  the  crust  is  a 
heterogeneous  body  whose  numerous  local  discontinuities 
comprise  the  thick  transition  zones  between  the  main  layers, 
in  particular,  between  the  crust  and  the  mantle. 

For  North  America,  only  isolated  instances  of  deep 
crustal  reflections  have  been  recorded  for  the  period  up 
to  1961.  Narans  et  al .  (1961)  performed  seismic  surveys 
in  Utah  in  an  effort  to  obtain  normal  incidence  reflections 
from  within  the  lower  crust.  In  order  to  distinguish 
between  random  noise  and  possible  reflections,  the  seismo¬ 
grams  were  analyzed  statistically  as  in  Germany.  They 
obtained  peaks  in  their  histograms  corresponding  to  discon¬ 
tinuities  at  8.5  and  26.3  kilometers.  However,  their 
results  were  based  on  only  ten  records. 

Dix  (1965)  applied  standard  reflection  prospecting 
methods  to  the  study  of  deep  crustal  structure  in  the 
Mojave  Desert.  To  eliminate  the  problems  of  multiple 
reflections  from  sedimentary  layers,  he  recorded  in  an 
area  of  thin  sedimentary  cover.  From  the  22  records,  he 
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found  strong  events  that  had  arrival  times  consistent  with 
that  of  the  Mohorovicic  discontinuity. 

The  first  deep  reflection  study  at  the  University 
of  Alberta  was  a  near  vertical  incidence  reflection  inves¬ 
tigation  in  southern  Alberta  (Kanasewich  and  Cumming,  1965). 
A  tapered  array  consisting  of  12  evenly  spaced  geophones 
over  134  meters,  with  an  additional  4  geophones  at  the 
central  location,  was  used.  The  distance  between  geophone 
units  was  293  meters.  This  type  of  array,  along  with 

multiple  holes,  was  found  to  be  an  effective  filter  for 

2 

attenuating  long  period  surface  waves.  From  a  T  versus 
2 

X  analysis,  the  average  velocity  to  a  strong  reflector  at 
11.4  seconds  was  found  to  be  6.37  kilometers  per  second. 

The  depth  calculated  to  this  reflector,  now  known  as  the 
Riel  discontinuity,  was  34  kilometers. 

Seismic  reflections  were  investigated  along  4  pro¬ 
files  comprising  some  90  kilometers  in  southern  Alberta 
(Clowes  et  al .  ,  1  968  ).  The  data  was  recorded  on  FM  analog 
magnetic  tape,  and  the  results  were  later  digitized. 

Through  power  spectra  calculations,  the  energy  of  the 
reflected  wavelets  was  found  to  be  in  the  5  to  15  hertz 
band.  Along  one  profile,  the  Riel  discontinuity  was 

located  at  11.6  seconds,  and  correlated  over  25  kilometers. 

2  2 

A  T  versus  X  analysis  yielded  an  average  velocity  of  6.2 
kilometers  per  second  to  a  depth  of  34  kilometers. 

Seismograms  obtained  from  the  studies  in  southern 
Alberta  were  later  investigated  to  determine  the  attenuation 


properties  of  the  crust,  and  the  characteristics  of  the 
reflecting  layers  (Clowes  and  Kanasewich,  1970).  The 
study  made  use  of  synthetic  seismograms  with  depth  and 
frequency  dependent  attenuations.  The  data  was  inverted 
by  comparing  the  autopower  spectra  of  particular  time 
intervals  on  the  field  records  with  the  spectra  computed 
from  the  synthetic  seismograms.  A  transition  zone  model 
comprised  of  alternating  high  and  low  velocity  sills 
satisfied  the  recorded  observations. 

During  the  summer  of  1969,  two  linear  arrays  were 
operated  through  a  joint  project  conducted  by  Princeton 
and  the  University  of  Alberta  (Perkins  and  Phinney,  1971). 
These  arrays  were  located  near  Riverton,  Wyoming.  The 
shot  and  geophone  patterns  were  as  described  by  Kanasewich 
and  Cumming  (1965).  Charge  weights  were  between  9  and  13 
kilograms,  and  were  detonated  at  30  meter  depths.  Data 
was  collected  on  FM  magnetic  tape  recorders  as  well  as  on 
a  galvanometer  camera.  Following  recording,  the  data  was 
digitized  so  that  signal  enhancement  techniques  could  be 
applied.  From  the  results,  good  reflections  were  noted  at 
11  and  35  ki 1 ometers  . 

A  near  vertical  incidence  reflection  survey  was 
conducted  in  north-central  British  Columbia  (Mair  and 
Lyons,  1976).  The  five  shot  points  were  separated  by  1.5 
kilometers,  and  the  42  recording  sites  were  spaced  at  250 
meters.  They  also  used  the  16  geophone  tapered  array  of 
Kanasewich  and  Cumming  (1965),  and  recorded  a  subsurface 


coverage  of  500  percent.  Two  7  channel  FM  tape  recorders 
were  used  to  record  the  data,  which  was  later  digitized 
with  a  sampling  interval  of  10  milliseconds.  Coherent 
energy  was  seen  to  exist  at  8  and  11  seconds.  The  event 
at  11  seconds  is  continuous  over  much  of  the  profile,  and 
is  believed  to  be  from  the  Mohorovicic  discontinuity.  All 
attempts  to  enhance  the  signal  from  the  Moho  by  the  common 
depth  point  (C.D.P.)  stack  met  with  failure,  and  they  sug¬ 
gest  that  this  is  the  result  of  the  crust-mantle  transition 
zone  being  a  many  layered  complex  of  high  and  low  velocity 
beds,  with  a  total  thickness  of  only  a  few  kilometers. 

In  1973,  an  additional  40  kilometers  of  reflection 
records  were  obtained  in  southern  Alberta  to  complement 
those  already  gathered  by  the  University  of  Alberta 
(Cumming  and  Chandra,  1975).  The  investigations  have  used 
power  spectra  analysis  and  predictive  deconvolution  of  the 
primary  discontinuities  within  the  crust  to  increase  our 
knowledge  of  their  nature. 

Several  reflection  profiles  were  shot  just  north  of 
Edmonton,  and  the  results  from  the  study  indicate  dips  within 
the  lower  crust  of  15  to  20  degrees  southeast  (Ganley  and 
Cumming,  1974).  Reflections  were  observed  at  20,  32  (Riel 
discontinuity)  and  35.5  (Mohorovicic  discontinuity)  kilo¬ 
meters  . 

A  near  vertical  incidence  reflection  profile  was 
carried  out  over  three  lines  in  Hardeman  County,  Texas 
(Oliver  et  al.,1976).  The  profile,  containing  368  kilo- 


meters  of  subsurface  coverage,  was  recorded  by  the  Vibro- 
seis  technique  utilizing  large  vibrators  and  low  frequency 
waves.  Deep  reflections  were  not  continuous  over  more 
than  a  few  kilometers.  Within  the  crust,  zones  of  low 
reflector  density  were  thought  to  be  the  result  of  plutons, 
and  curvature  of  the  reflections  to  indicate  deep  folded 
structures.  The  Mohorovicic  discontinuity  was  located  at 
14  to  15.4  seconds. 

Smithson  and  Brown  (1977)  propose  a  model  for  the 
lower  crust  based  on  geological,  geochemical  and  geophysi¬ 
cal  information.  Their  geophysical  data  comes  from  both 
refraction  and  reflection  studies.  Discussion  of  the  field 
procedures  and  the  initial  results  of  the  refraction  study 
are  reported  by  Smithson  and  Shire  (1975).  Smithson 
and  Brown  suggest  that  the  lower  crust  is  a  heterogeneous 
body,  with  lateral  heterogeneity  being  shown  by  the  seismic 
wave  anomalies  across  the  arrays.  The  number  of  reflectors 
in  the  lower  crust,  and  the  presence  of  composite  reflec¬ 
tions,  are  best  interpreted  as  layering  found  in  metamor- 
phic  rocks.  They  suggest  that  the  lower  crust  consists  of 
granite  and  syenite  gneiss,  as  well  as  amphibolite,  which 
are  interlayered  and  intruded  by  gabbro.  The  characteris¬ 
tic  discontinuity  of  deep  reflections  are  thus  the  result 
of  disruptions  of  the  layers,  tight  folding,  changes  in 
dip  and  layer  thickness,  and  intrusions. 
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1.3.  The  Reflection  Method 

The  seismic  reflection  technique  was  developed  for 
detailing  the  structure  of  the  Earth  by  measuring  the 
times  necessary  for  a  seismic  pulse  to  travel  from  a 
source  to  a  reflector,  and  back  to  the  surface  again.  In 
seismic  operations,  a  signal  is  induced  into  the  Earth 
either  by  the  detonation  of  a  charge  in  a  shot  hole,  or  by 
some  mechanical  means  such  as  Vibroseis  or  Di  nose  is 
(Dobrin,  1976)  (figure  1.2).  Upon  encountering  various 
discontinuities  in  the  Earth,  some  of  the  propagating 
energy  is  reflected  back  to  the  Earth's  surface.  The 
quantity  of  energy  that  is  reflected  depends  primarily  on 
the  velocity  and  density  distribution  found  in  the  Earth. 
This  reflected  energy  is  collected  along  arrays  of  geo¬ 
phone  groups,  each  group  consisting  of  several  geophones. 
From  the  geophones,  the  reflected  signals  are  time  multi¬ 
plexed,  amplified,  converted  to  digital  form,  and  written 
onto  magnetic  tape  (Mercado,  1973). 

Through  the  use  of  digital  computers,  the  reflected 
signals  are  brought  into  proper  time  alignment,  so  as  to 
simulate  normal  incidence  seismic  reflection  profiling. 
This  data  is  then  passed  through  various  filters  to 
enhance  the  signal  relative  to  the  noise.  The  subject  of 
this  thesis  is  first  to  discuss  those  techniques  used  to 
bring  reflected  events  into  proper  time  alignment,  and 
then  to  investigate  which  methods  best  enhance  the  signal 
with  respect  to  the  noise. 
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Figure  1.2:  The  reflection  technique. 


CHAPTER  2 
THE  CRUSTAL  STUDY 


2.1.  Introduction 

The  transition  zone  between  the  Churchill  (figure 
2.1)  and  Superior  structural  provinces  separates  rocks 
which  differ  not  only  in  age,  but  in  lithology  and  struc¬ 
tural  trends  as  well.  In  northeast  Manitoba,  structural 
trends  of  the  Superior  Province  are  expressed  in  the  east- 
west  orientations  of  metavolcanic  and  metasedimentary  belts 
(Goodwin,  1972).  These  supracrustal  rocks  are  Archean  in 
age  ( K -  A  r  radiometric  ages  of  about  2500  million  years), 
and  lie  in  synformal  bands  which  separate  the  much  more 
extensive  granites  and  gneisses.  The  volcanic  rich  belts 
(greenstone  belts)  consist  of  varied  volcanic  assemblages 
with  minor  amounts  of  sediments,  and  have  been  extensively 
intruded  by  granodi ori ti c  plutons.  As  one  moves  away  from 
the  bounding  granitic  batholiths,  the  metamorphic  grade 
decreases  from  amphibolite  facies  to  greenschist  facies. 

On  the  other  hand,  the  metasedimentary  belts  contain 
assorted  metasediments  along  with  schist,  migmatite,  and 
paragneiss.  The  metamorphic  grade  is  fairly  uniform  at 
the  amphibolite  grade,  but  greenschist  and  granulite 
facies  may  be  found  . 

Within  the  Manitoba  portion  of  the  Superior  province, 
strong  east-west  magnetic  and  gravitational  trends  exist. 
Positive  Bouguer  anomalies  occur  over  the  supracrustal 
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belts  and  granulite  rocks,  with  negative  anomalies  over- 
lying  the  granitic  rocks.  The  stronger,  linear  magnetic 
anomalies  coincide  with  the  greenstone  belts  and  granulites, 
with  the  low  intensity  areas  corres pondi ng  to  the  para- 
gneisses  and  granites. 

The  Churchill  province  in  northern  Saskatchewan  is 
underlain  mainly  by  granitic  and  gneiss ic  rocks.  Curved 
belts  of  metamorphosed  and  deformed  lower  Proterozoic 
(K-Ar  radiometric  ages  around  1800  million  years)  supra- 
crustal  rocks  are  found.  These  supracrustal  rocks  are 
mainly  metasediments,  migmatites  and  slightly  metamorphosed 
volcanic  rocks.  The  predominant  structural  trend  in  the 
western  half  of  the  province  is  northeast,  but  it  slowly 
changes  to  the  southeast  as  one  moves  towards  the  east. 

Both  the  magnetic  and  gravity  anomalies  found  within  the 
Saskatchewan  portion  of  the  Churchill  province  are  highly 
erratic,  and  do  not  show  any  strong  trends. 

Both  geological  and  geophysical  information  have 
been  used  to  define,  and  investigate  the  boundary  between 
the  Churchill  and  Superior  structural  provinces  (Lee, 

1977).  Early  attempts  at  defining  the  boundary  are  found 
in  Bell  (1971a).  Burwash  et  al .  (  1  962  )  first  proposed  the 
boundary  to  be  a  transition  zone,  which  was  later  affirmed 
by  Bell  (1971b).  In  his  paper,  Bell  defined  the  transition 
zone  in  northern  Manitoba  as  the  Pikwitonei  subprovince. 
However,  Rb-Sr  radiometric  dating  (Cranstone  and  Turek, 

1976)  in  the  boundary  zone,  along  with  evidence  of 


Pikwitonei  rocks  being  traced  into  the  Wabowden  subprovince 
(Green  et  al .  ,  1977),  suggest  that  this  is  insufficient. 
Thus,  the  Wabowden  subprovince  is  now  included  in  the 
boundary  zone.  In  this  boundary  zone,  rocks  of  the  Pik¬ 
witonei  subprovince  are  Archean  in  age  (Cranstone  and 
Turek,  1976),  while  in  the  Wabowden  subprovince,  ages  are 
typically  lower  Proterozoic.  Rocks  of  the  Pikwitonei 
subprovince  are  predominantly  granulite  gneisses  and 
charnockites  (contain  amounts  of  hypersthene)  (Bell, 

1971b).  Migmatitic  gneisses  predominate  in  the  Wabowden 
subprovince.  The  western  edge  of  this  subprovince  comprises 
the  Thompson  Nickel  Belt  (Lee,  1977). 

In  the  north,  the  boundary  between  the  Wabowden  and 
Pikwitonei  subprovinces  is  a  major  fault  zone  (  Assean  Lake 
fault  zone),  as  described  by  Bell  (1971b).  Cranstone  and 
Turek  (1976)  suggest  that  the  boundary  zone  is  a  meta- 
morphic  transition  zone  under  the  Phanerozoic  sediments. 

Three  lines  of  evidence  indicate  the  extension  of 
the  Church i 1 1 -Superi or  boundary  zone  under  the  Phanerozoic 
sediments  of  the  Interior  Plains.  Geochronological  data 
from  core  samples  show  that  both  the  Churchill  and 
Superior  provinces  extend  into  southern  Manitoba, 
Saskatchewan  (Burwash  et  al .  ,  1962  ;  Goldich  et  al .  ,  1966), 
and  into  North  Dakota. 

The  Bouguer  gravity  map  of  eastern  Saskatchewan  and 
western  Manitoba  (figure  2.2)  shows  that  the  Nelson  River 


Figure  2.2:  Gravity  map  of  southern  Saskatchewan  and  southern 

Manitoba  with  R=Regina  and  W=Winnipeg  (after  Lee, 1977). 
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Figure  2.3:  Magnetic  map  of  southern  Saskatchewan  (after  Lee, 1977) 
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gravity  high  can  be  followed  through  the  northern  parts  of 
Manitoba  into  southeast  Saskatchewan.  In  northern  Manitoba, 
Gibb  (1968a,  1968b)  has  shown  that  this  gravity  high  is 
associated  with  the  high  density  granulites  of  the 
Pikwitonei  subprovince.  On  the  other  hand,  the  Wabowden 
subprovince  is  associated  with  a  gravity  low. 

The  boundary  zone  in  northern  Manitoba  is  also  well 
expressed  in  aeromagnetic  maps  (figure  2.3).  Granulite 
gneisses  of  the  Pikwitonei  subprovince  are  expressed  in  a 
configuration  of  oval  shaped  high  and  low  magnetic  anomalies 
(Bell,  1971b;  Kornik  and  Maclaren,  1966).  In  comparison 
with  this,  the  migmatitic  gneisses  of  the  Wabowden 
subprovince  are  exemplified  by  high  density,  elongated 
anomalies,  which  trend  to  the  northeast.  Green  et  al . 

(1977)  indicate  that  these  elongated  anomalies  may  be 
recognized  in  southeast  Saskatchewan  and  southwest  Manitoba. 
On  this  basis,  Green  has  delineated  the  boundary  zone 
under  the  Phanerozoic  sediments  to  the  Canada-United  States 
border. 

The  overall  picture  of  the  Churchill-Superior 
transition  zone  as  suggested  by  geophysical  and  geological 
data  is  that  this  zone  may  be  an  ancient  Precambrian 
suture  zone  (  Dewey  and  Burke,  1973).  This  boundary  zone 
is  associated  with  two  important  geological  deposits 
(Green  et  al .  ,  1  97  7  ).  In  northern  Manitoba,  nickel 
is  mined  from  the  Thompson  Nickel  Belt,  one  of  the 


23 


world's  largest  nickel  deposits.  In  the  southern  regions 
of  Saskatchewan  and  Manitoba,  the  Phanerozoic  sediments  over- 
lying  the  transition  zone  have  entrapped  large  quantities  of 
oil  and  gas.  Though  much  of  the  cratonic  movement  between 
the  Churchill  and  Superior  provinces  would  have  ceased  by 
the  Phanerozoic,  these  traps  may  be  the  result  of  slight 
readjustments  between  the  two  plates,  due  to  persisting 
stresses,  resulting  in  deformation  of  the  overlying  sediments. 
Thus,  locating  these  oil  reserves,  and  perhaps  other  nickel 
deposits,  would  be  far  easier  if  the  geological  evolution  of 
the  rocks  within  the  boundary  zone  were  known.  Therefore, 
it  was  decided  that  a  detailed  geophysical  study  of  the 
Churchill-Superior  boundary  zone,  including  seismic,  gravity 
and  magnetic  studies,  would  be  undertaken  to  increase  our 
knowledge  of  this  boundary  zone. 

2.2.  The  Project 

During  1975  and  1976,  at  the  suggestion  of  Dr.  E.R. 
Kanasewich,  a  loose  association  of  seismologists  from  the 
universities  of  British  Columbia,  Alberta,  Manitoba, 

Toronto  and  Western  Ontario  and  the  Earth  Physics  Branch 
of  the  Department  of  Energy,  Mines  and  Resouces  met 
several  times  to  form  a  Canadian  Crustal  Studies  Group 
(CCSG).  In  July  1977,  the  CCSG  performed  its  first  seismic 
reflection  survey  in  southern  Saskatchewan  and  Manitoba 
under  the  direction  of  Dr.  A.G.  Green.  The  project  con¬ 
sisted  of  a  wide  angle  reflection/refraction  profile, 
along  with  a  near  vertical  incident  reflection  survey. 
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The  main  purpose  of  this  survey  was  to  increase  our  know¬ 
ledge  of  the  transition  zone  between  the  Churchill  and 
Superior  structural  provinces. 

A  number  of  factors  influenced  the  decision  to  have 
the  survey  in  southern  Manitoba  as  opposed  to  its  northern 
regions,  where  the  boundary  zone  is  visible.  First,  there 
was  evidence  for  both  deep  and  shallow  reflecting  layers 
in  the  crust,  with  good  control  on  depths  and  nature  of 
the  shallow  layers  from  oil  company  drill  holes.  The  topo¬ 
graphy  was  extremely  flat,  with  a  highly  unconsolidated 
surface  layer  which  allowed  for  easy  drilling  of  shot  holes. 
Also,  there  was  a  complete  network  of  roads  available  which 
yielded  easy  access  for  drilling,  shooting  and  recording. 

The  combined  survey  consisted  of  a  245  kilometer 
north-south  and  a  242  kilometer  east-west  refraction  sur¬ 
vey,  and  a  84.54  kilometer  east-west  near  vertical  incident 
reflection  survey  (figure  2.4).  The  reflection  survey  was 
conducted  over  three  lines  which  were  44.17,  19.31  and 
21.06  kilometers  in  length  (figure  2.5),  with  an  average 
subsurface  coverage  of  400  percent  (figure  2.6).  Alto¬ 
gether,  there  were  58  shots  detonated  at  an  average  depth 
of  18  meters  in  holes  that  were  15  centimeters  in  diameter. 
The  explosive  used  consisted  of  60  percent  geogel  in  46  by 
13  centimeter  sticks.  For  the  spread,  there  were  a  total 
of  47  takeouts.  With  a  takeout  distance  of  292.55  meters, 
this  yielded  a  spread  length  of  12.87  kilometers  (figure 
2.7).  Each  trace  was  recorded  with  a  linear  array  of  nine 
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Figure  2.4:  Map  of  the  project. 
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Figure  2.5:  The  reflection  survey. 
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Figure  2.6:  Histograms  of  subsurface  coverage. 


L-15A,  600  ohm-10  hertz  geophones  on  a  73  meter  A30  DDC 
stringer.  At  the  shot  hole,  an  uphole  geophone  was  also 
recorded  . 

The  shooting  sequence  is  depicted  in  figure  2.8.  In 
most  cases,  the  shot  holes  were  located  at  every  fifth  or 
sixth  receiver.  After  four  shots,  the  complete  spread  was 
moved  westward  6.44  kilometers.  For  every  shot  the  detona 
tion  circuit  was  connected  to  the  recording  cable  of  the 
Manitoba  system,  and  to  a  chronometer  unit  that  detonated 
the  shot  (Green  et  at,,  1977).  The  WWVB  radio  signal  was 
synchronized  with  the  chronometer  unit  so  as  to  provide 
accurate  timing  of  the  shot  instant.  The  chronometer  unit 
was  also  connected  to  the  time  break  cables  which  recorded 
a  tone  break  when  the  shot  was  detonated.  This  tone  break 
was  transmitted  by  radio  to  the  other  systems.  The  shot 
instant  was  taken  to  be  the  time  of  the  tone  break,  with 
the  delay  in  detonation  from  the  0.0  second  mark  on  the 
WWVB  signal  and  the  tone  break  being  about  42  milliseconds 

In  recording  the  data,  three  trucks  were  used, 
corresponding  to  the  three  systems.  The  University  of 
Manitoba  employed  a  DS1590  automatic  gain  varying  digital 
system,  and  recorded  at  a  sample  rate  of  one  millisecond. 
They  recorded  traces  11  through  34  along  the  spread.  The 
University  of  Saskatchewan  had  a  programmed  gain  analog 
recorder.  Traces  34  to  45  were  recorded  with  their  system 
A  fixed  gain  digital  system  which  sampled  at  5.6  milli¬ 
seconds  was  employed  by  the  University  of  Alberta.  This 
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Figure  2.8:  The  shooting  sequence  (modified  from  Green  et  al 


system  was  responsible  for  recording  traces  1  to  1 1  .  Trace 
1  was  recorded  at  the  most  easterly  receiver,  and  the  two 
overlapping  traces  allowed  a  tie  between  the  three  systems. 

Since  much  preliminary  work  had  to  be  done  on  the 
data  before  it  could  be  put  into  a  form  suitable  for 
analysis,  it  is  worthwhile,  at  this  time,  to  discuss  our 
system,  and  its  mode  of  data  collection. 

2.3.  The  Digital  Recording  System 

The  University  of  Alberta's  recording  system  is  a 
digital  recording  system.  With  the  system,  14  channels  of 
information  may  be  digitally  recorded  onto  an  800-bytes 
per  inch,  nine-track  synchronous  tape  recorder.  The  tape 
velocity  is  6.25  inches  per  second,  yielding  a  5000  bytes 
per  second  transfer  rate.  During  the  survey,  channel  7 
was  used  to  record  the  WWVB  radio  signal,  channel  14  to 
record  the  tone  break,  and  channels  1  through  5  and  8 
through  13  were  used  to  record  the  seismic  signals.  These 
channels,  excluding  the  dead  channel  6,  were  also  recorded 
photographically. 

Figure  2.9  details  the  complete  recording  system. 

The  tape  recorder  is  a  model  7830-9  write-only  transport, 
and  the  multiplexer,  A  to  D  converter  and  5  kilohertz 
crystal  clock  are  contained  within  a  modified  model  120 
Data  Acquisition  System  (All  sop  et  al.>  1972).  One  to 
twenty  channels  may  be  recorded  with  a  maximum  of  4096 
samples  per  channel.  Several  records  may  be  sequentially 
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Figure  2.9:  The  digital  recording  system(after  Allsop  et  al  ,1972) . 
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written,  with  an  interrecord  gap  of  0.096  seconds. 

The  signal  is  converted  from  analog  to  digital  form 
by  a  14  bit  converter  (13  bits  plus  sign)  with  a  dynamic 
range  of  84  decibels.  Since  two  8-bit  bytes  are  required 
to  store  each  sample,  2500  samples  are  recorded  per 
second.  During  the  survey  we  recorded  14  channels,  and 
thus  the  sampling  rate  is  14/2500  (or  0.0056)  seconds. 

This  yields  a  Nyquist  frequency  of  89.3  hertz.  The  length 
of  the  record  written  was  4096  samples  in  integer  *2 
words,  and  thus  we  recorded  22.94  seconds  of  data.  With 
the  collection  of  14  channels,  that  results  in  a  block  length 
of  114,688  bytes. 

The  amplifier  circuit  is  presented  in  figure  2.10, 
and  is  discussed  in  some  detail  by  Ganley  (1973).  Briefly, 
though,  it  may  be  described  as  follows.  A  fixed  gain  40  db 
preamplifier  is  driven  by  an  18  db  gain  input  transformer. 
This  preamplifier  feeds  the  main  amplifier  through  a  two- 
pole  high-pass  Bessel  filter.  The  Bessel  filter  yields 
zero-phase  shift,  and  may  be  bypassed.  Output  from  the 
main  amplifier  (which  has  switchable  gain,  with  increments 
of  10  db  from  0  db  through  60  db)  is  passed  through  a  four- 
pole  Bessel  aliasing  filter  with  a  50  hertz  corner.  The 
gain  and  phase  response  of  the  amplifier  are  discussed 
fairly  completely  by  Dick  ins  (  1  973  ). 

Figure  2.11  exhibits  one  of  the  photographically 
recorded  seismic  sections.  Again,  channel  7  is  the  WWVB 
signal  and  channel  14  is  the  tone  break.  Due  to  the  high 
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.10:  The  amplifier  circuit  (after  Allsop  et  al 


972) 


gain  setting  necessary  for  recording  deep  reflections, 
the  information  in  the  early  section  is  clipped.  This 
same  section,  as  output  from  the  digital  recording  system, 
is  seen  in  figure  2.12  (the  WWVB  radio  signal,  the  dead 
channel  6,  and  the  tone  break  are  not  shown).  The  separa¬ 
tion  of  the  individual  channels  in  this  display  makes 
visual  interpretation  of  the  seismic  section  much  simpler. 
Also,  it  is  possible  to  compress  the  channels  into  a  much 
smaller  length  as  compared  to  the  photographic  record  which, 
on  the  average,  is  about  three  meters  in  length. 
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Figure  2.11:  Photographic  record  for  shot  123. 
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Figure  2.12:  Digital  output  plot  for  shot  123(after  redigitization). 


CHAPTER  3 
DIGITAL  ANALYSIS 


3.1.  Preliminary  Analysis 

The  reflection  seismic  profile  was  gathered  on  three 
different  recording  systems  in  the  common  depth  point 
(C.D.P.)  format.  To  be  able  to  properly  interpret  the 
data,  all  of  it  must  be  contained  on  one  tape,  and  in  an 
order  suitable  for  C.D.P.  studies.  Thus,  our  data  had  to 
be  first  made  compatible  with  that  of  the  other  universi¬ 
ties. 

The  first  step  in  the  treatment  of  the  field  recorded 
data  was  to  demultiplex  it,  and  write  it  onto  a  single 
tape.  Each  file  of  this  first  tape  consisted  of  the  four¬ 
teen  channels,  each  containing  4096  samples,  recorded  for 
each  shot.  Next,  the  data  was  converted  from  integer  *2 
words  to  real  *4  words,  and  then  redigitized  so  that  all  of 
the  universities'  information  had  the  same  sampling  inter¬ 
val.  The  sampling  interval  chosen  was  four  milliseconds, 
yielding  seismic  traces  with  250  samples  per  second.  In 
the  course  of  his  research,  D.  Ganley  has  written  a  pro¬ 
gram  that  will  resample  a  time  series  from  one  sampling 
interval  to  another.  His  algorithm  makes  use  of  the  spline 
interpolation  theory  found  in  Pennington  (1970).  Basically, 
the  program  connects  each  pair  of  adjacent  points  with  a 
section  of  a  third  degree  polynomial,  such  that  the  first 
and  second  derivatives  at  each  point  are  continuous. 
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Also,  at  each  end,  the  second  derivative  is  a  linear  extra¬ 
polation  of  the  value  at  (  X  (  2  )  ,  Y ( 2  )  )  and  (X(N-1  )  ,Y(N-1  )) 
that  is,  the  third  derivative  is  continuous  at  these 
points  (N  is  the  number  of  points  in  the  original  time 
series;  where  X  and  Y  are  coordinates  of  points  in  the  input 
time  series).  The  result  of  this  redigitization  may  be 
seen  in  figure  3.1.  This  illustration,  in  the  upper  half 
(traces  a  and  b),  compares  the  input  and  output  traces; 
the  output  series  now  containing  5734  points.  In  the 
clipped  portions  of  the  output  trace  (b),  we  notice  that 
small  spikes  have  been  generated  on  the  edges  of  the 
blocks.  These  peaks  are  very  high  frequency,  and  are 
eliminated  through  frequency  filtering  (trace  c  of  figure 
3.1).  Figure  3.2  shows  the  output  power  of  the  trace 
before  (a)  and  after  (b)  redigitization.  The  plots  are 
quite  similar,  except  for  a  large  peak  at  118  hertz  for 
the  redigitized  trace.  This  is  due  to  the  small  spikes 
generated  in  the  redigitization  process.  We  should  also 
note,  at  this  time,  that  the  new  Nyquist  frequency  is  125 
hertz  . 

Following  redigitization,  it  was  necessary  to  locate 
all  sensors  that  were  connected  with  reverse  polarity  in 
the  field.  This  required  scanning  the  analog  monitor 
records  of  all  the  shots  and  picking  out  traces  with  first 
motion  opposite  to  the  other  traces  on  the  record.  Only 
a  few  sensors  were  reversed.  At  this  time,  it  was  conven¬ 
ient  to  make  sure  that  all  of  the  records  had  the  same 
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Figure  3.1:  Comparison  of  a  seismic  trace  before  and 

after  redigitization.  Trace  (a)  is  the 
original  seismic  channel.  Trace  (b)  is  the 
same  channel  after  redigitization  and  the 
addition  of  450  points  of  noise.  It  has 
also  been  divided  through  by  four.  This 
channel,  after  8  to  40  hertz  band-pass 
filtering,  is  seen  as  trace  (c). 
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Figure  3.2:  Comparison  of  the  log  of  the  power  for  the  seismic 

channel  of  figure  3.1  before  and  after  redigitization. 
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block  length  and  that  the  shot  instant  was  superimposed  on 
all  traces.  It  was  decided  that  each  trace  would  contain 
450  points  of  noise  before  the  shot  instant  and  5100 
points  of  data,  yielding  a  block  containing  5550  samples. 

To  accurately  pick  the  shot  instant  for  each  trace,  a 
relatively  simple  procedure  was  adopted.  Each  file  of 
the  tape,  containing  the  fourteen  channels  for  each  shot, 
were  plotted  out  and  compared  with  the  corresponding  field 
photographic  recording.  Using  the  photographic  record, 
the  shot  instant  was  picked  exactly,  knowing  that  the 
blaster  introduced  a  42  millisecond  delay  after  the  0.0 
time,  which  was  shown  on  the  WWVB  radio  signal  immediately 
below.  Then,  knowing  exactly  when  the  shot  instant 
occurred,  times  between  various  points  along  the  traces 
and  the  shot  tone  were  measured  on  the  paper  record.  All 
of  these  time  intervals  were  accurate  to  about  one  milli¬ 
second.  These  same  time  intervals  were  then  measured  on 
the  plot,  and  in  most  cases,  were  within  four  milliseconds 
of  the  times  measured  on  the  paper  record.  However,  in 
some  cases,  the  tone  break  was  not  recorded  digitally  very 
well,  and  so  it  was  necessary  to  use  the  times  measured  on 
the  paper  record  to  extrapolate  the  point  corresponding  to 
the  shot  instant  on  the  plot.  It  was  possible,  in  all 
cases,  to  pick  the  shot  instant  to  an  accuracy  of  one 
point. 

The  next  step  in  the  preparation  of  the  tape  was  to 
filter  the  data.  Since  some  of  the  cables  were  near 
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power  lines  (as  evidenced  by  the  peak  in  the  power  spectrum 
at  60  hertz  in  figure  3.2),  an  upper  cutoff  frequency  of 
55  hertz  was  chosen.  A  lower  cutoff  frequency  of  3  hertz 
was  also  chosen.  Digital  filtering  of  the  data  will  be 
discussed  in  more  detail  in  section  3.4,  but  basically 
the  data  were  passed  through  a  zero-phase  shift,  eight-pole 
recursive  algorithm  based  on  a  Butterworth  function  with 
attenuation  outside  the  passband  being  96  db  per  octave. 

The  data  were  also  filtered  with  an  8  to  40  hertz  filter. 

Finally,  we  had  to  convert  the  data  back  into 
integer  *2  words.  The  final  tape  consisted  of  nine  files; 
three  files  containing  the  three  lines  of  raw  data,  three 
files  with  the  3  to  55  hertz  filtered  data,  and  three  files 
with  the  8  to  40  hertz  filtered  data.  The  files  contained 
only  the  data  channels  for  each  shot,  and  each  trace  had 
a  label  containing  information  on  the  shot  and  receiver 
locations.  Since  each  trace  made  up  a  block,  the  block 
length  was  11,000  +  100  (label)  bytes.  This  tape  was  sent 
to  Manitoba  where  the  data  from  all  three  of  the  universi¬ 
ties  was  sorted  according  to  subsurface  coverage  point, 
and  compiled  onto  one  tape.  This  new  tape  contains  three 
files  corresponding  to  the  three  lines,  with  the  informa¬ 
tion  being  stored  in  integer  *2  words.  The  Manitoba  and 
Alberta  traces  have  had  a  3  to  55  hertz  filter  applied  to 
them,  and  the  Saskatchewan  traces  have  been  subjected  to 
a  0  to  55  hertz  filter.  For  the  sake  of  safety,  a  copy 
of  this  tape  was  made. 
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3.2. 


Static  Corrections 


To  process  C.D.P.  data  properly,  all  of  the  traces 
corresponding  to  a  particular  subsurface  point  must  be 
brought  into  proper  time  alignment.  Timing  corrections 
are  of  two  forms.  Static  corrections  are  based  on  spatial 
variations  of  velocity  and  on  the  variations  of  elevation, 
and  dynamic  corrections  are  time  dependent  due  to  varia¬ 
tions  in  ray  parameter,  and  will  be  discussed  in  section 
3.5.  Each  static  correction  is  the  sum  of  an  elevation 
correction  plus  a  near  surface  low  velocity  (weathered)  cor¬ 
rection  (Telford  et  al . ,  1976).  Because  of  the  anomalously 
low  velocity  in  this  layer,  any  changes  in  layer  thickness 
will  overly  affect  the  travel  times  of  the  reflected  waves, 
and  thus,  it  must  be  corrected  for  (Dobrin,  1976).  The 
elevation  correction  placed  all  shot  points  and  receivers 
at  a  common  datum  of  380  meters  above  sea  level.  By 
measuring  the  time  interval  between  the  tone  breaks  and 
the  first  breaks  on  the  uphole  geophones  at  all  of  the 
shot  points,  the  velocity  of  the  seismic  waves  in  the 
weathered  layer  was  determined  to  be  1.89  kilometers  per 
second.  For  those  receiver  locations  that  were  not  used 
as  shot  points,  the  weathering  corrections  were  interpol¬ 
ated.  These  corrections  were  applied  to  the  entire  tape, 
and  in  conjunction  with  this,  the  data  was  converted  into 
real  *4  format . 

My  thesis  is  primarily  concerned  with  investigating 
the  various  techniques  used  to  interpret  seismic  reflection 
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data,  and  so  only  small  sections  of  the  tape  were  worked 
with  at  any  one  time.  With  the  data  being  stored  in 
C.D.P.  format,  each  subsurface  point  could  comprise  data 
from  one,  two  or  all  three  of  the  universities.  The 
sections  of  data  that  I  chose  involved  traces  recorded  by 
both  the  universities  of  Alberta  and  Manitoba.  Since  each 
university  independently  picked  the  shot  instants  for  each 
shot,  they  were  probably  picked  differently.  However,  in 
order  to  obtain  proper  time  alignment  between  traces,  the 
shot  instant  must  occur  at  the  same  time  for  each  trace 
corresponding  to  each  subsurface  point.  This  was  easily 
assured  by  comparing  the  two  overlapping  traces  for  each 
shot  (traces  11  and  12).  A  program  was  written  which  cross 
correlates  two  seismic  traces,  and  then  plots  them  both 
out.  By  visually  comparing  the  two  traces,  and  utilizing 
the  cross  correlation  coefficients,  I  was  able  to  give 
the  traces  recorded  by  the  two  universities  the  same  shot 
instant.  An  example  of  this  procedure  is  shown  in  figures 
3.3  and  3.4.  The  first  figure  exhibits  the  two  overlapping 
traces  for  shot  34A.  The  lower  trace  was  recorded  by  the 
Univerisity  of  Alberta,  and  was  used  as  the  standard  trace 
in  the  cross  correlation  program.  Figure  3.4  is  the  cross 
correlation  curve  showing  that  trace  11  is  ahead  of  trace 
12,  that  is,  the  Manitoba  data  is  two  lag  points  behind 
the  Alberta  data  for  shot  34A.  In  all  cases,  I  assumed 
that  I  had  picked  the  shot  instances  correctly,  and  thus, 
Manitoba's  traces  were  adjusted  accordingly.  Generally 
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Figure  3.4:  Cross  correlation  curve  between  traces  11  and  12  for 
shot  34A.  Trace  11  is  the  reference  channel. 


though,  the  shot  instances  were  picked  within  three  points 
of  one  another . 


3.3.  Power  Spectral  Analysis 

In  order  to  determine  the  frequency  content  of  the 
reflection  seismograms,  a  computer  program  was  written 
which  computed  power  spectrums.  The  program  first  selects 
a  particular  time  window  of  a  trace  to  analyze,  and  then 
stores  it  in  a  vector.  To  prevent  power  from  being  gener¬ 
ated  at  all  frequencies  by  the  presence  of  sharp  end  dis¬ 
continuities,  the  first  and  last  ten  percent  of  the  data 
points  are  tapered  to  the  mean  value,  which  is  often  zero. 
The  taper  consists  of  a  pair  of  cosine  bells  (Kanasewich, 
1975).  For  N  data  points,  the  nth  point  is  multiplied  by 

Wn  =  Tj-  (1  +  cos  ---  —Ll)  -(L+M)  <  n  <  -  L 


=  (1  +  cos  11  ( nH~ 1 }  )  L  <  n  <  L  +  M 

The  period  of  the  cosine  bell  is  determined  by  M. 

Following  this,  the  digitized  seismic  trace  is  trans¬ 
formed  into  the  frequency  domain  through  the  use  of  the 
discrete  Fourier  transform,  which  may  be  represented  as 


‘ 


F(W)  =  I’'  f(n)e-27TiWn/N 
n  =0 


where 


W=0 , . . . ,  N- 1 

3.3.2 


f(n)  =  nth  digital  point 

The  algorithm  used  to  calculate  the  Fourier  transform  is 

based  on  the  fast  Fourier  transform  (FFT)  algorithms  of 

Cooley  and  Tukey  (1965).  To  employ  the  algorithm,  the 

N  data  points  are  factored  into  K  products.  Then  the 

computation  of  K  shorter  Fourier  transforms  combined  with 

a  series  of  recursive  multiplications  of  these  transforms 

(Kanasewich,  1975),  will  yield  the  Fourier  transform  of  the 

seismic  trace.  The  number  of  operations  in  the  calculation 

2 

reduces  from  N  to  N  Log2N.  In  the  algorithm  used,  N, 
the  number  of  data  points,  is  factored  into  products  of  two, 
and  thus,  the  number  of  data  points  must  be  a  power  of  two. 
This  is  facilitated  by  padding  both  ends  of  the  tapered 
trace  with  zeros  to  the  nearest  power  of  two  (N1).  The 
addition  of  the  zeros  to  the  trace  does  not  alter  its  rela¬ 
tive  power  spectrum,  but  does  smooth  the  spectral  window. 

The  algorithm  used  to  obtain  the  fast  Fourier  transform  of 
the  data  was  found  in  Claerbout  (1976).  For  a  complete 
discussion  of  Fourier  transforms,  including  the  FFT,  one 
is  referred  to  Kanasewich  (1975),  Claerbout  (1976), 

Dobrin  (1976),  and  to  Lindseth  (1967). 

A  Daniel!  spectral  estimate  is  now  computed  by  the 


program.  An  ideal  window  would  be  one  whose  frequency 
response  is  rectangular  (figure  3.5),  and  the  Daniell  win¬ 
dow  approaches  this  ideal  quite  well.  This  window  is 
computed  directly  on  the  periodogram  (Kanasewich,  1975)  by 


P(W)  = 


IT _ 

N  2m  +  1 


m 


1  l  F(W-j)F*(W-j)  W  =  0,1 . 


J  =-m 


3.3.3 


The  factor 


1 


2m  +  l 


normalizes  the  area  under  the 


window  to  unity,  and  N 1 / N  keeps  the  amplitude  of  the 
power  spectrum  constant  for  all  values  of  N'.  Since  the 
data  are  real  , 


f (n)  =  f*(n) 


3.3.4 


and  thus 


F(W)  =  F* ( - W ) 


3.3.5 


The  Fourier  transform  function  is  cyclic,  and  so 


F  (  - W  )  =  F(N-W) 


3.3.6 


Using  equations  3.3.3,  3.3.5  and  3.3.6,  the  Daniell 
spectral  estimate  may  be  written  as 


m 


P(W)  =  x  2iTT  £  F(W-j)(N-W+j)  W=0 


j=-m 


rr 

2 


9  •  •  •  9 


3.3.7 


m 


m 


W(go)=1  _T^4(j0  4  7T^ 

~  m  m 

W(oj)=0  otherwise 


Figure  3.5:  The  ideal  window  (after  Kanasewich,1975) . 
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The  power  of  various  seismic  channels  was  calculated 
over  the  time  window  from  6  to  18  seconds  to  determine 
their  frequency  content.  In  most  cases,  m  was  set  equal 
to  two.  Figure  3.6  exhibits  the  power  output  of  two  traces. 
From  the  plots,  we  can  see  that  the  main  power  lies  in  the 
10  to  30  hertz  range.  Assuming  a  crustal  velocity  of  6.5 
kilometers  per  second,  this  suggests  that  the  seismic  energy 
from  the  crust  and  mantle  is  propagating  with  wavelengths 
of  200  to  450  meters.  Thus,  the  resolution  obtainable  with 
this  data  will  be  of  the  order  of  200  to  450  meters. 

On  the  power  plots  of  raw  data,  peaks  are  seen  at 
about  3  and  60  hertz,  and  sometimes  at  even  higher  frequen¬ 
cies  (figure  3.2).  The  very  high  frequency  peaks  (with 
exclusion  of  the  118  hertz  peak)  are  likely  due  to  wind 
noise,  which  often  comes  in  with  a  frequency  of  about 
85  hertz.  The  60  hertz  peak  is  the  result  of  geophones 
being  located  near  power  lines,  and  the  3  hertz  peak  is 
likely  due  to  ground  roll.  Ground  roll  (roller)  is 
most  likely  the  result  of  surface  Rayleigh  waves,  and  it 
is  most  strongly  generated  where  the  near  surface  material 
is  highly  unconsolidated  and  of  very  low  velocity 
(Mateker,  1965).  However,  the  proper  arrangement  of 
geophones  and  the  use  of  multiple  holes  could  have 
severely  attenuated  the  ground  roll  (Kanasewich  and 
Cumming,  1965),  and  greatly  improved  the  quality  of  the 
records.  For  a  complete  discussion  on  the  theory  behind 
such  patterned  arrays,  along  with  their  amplitude  responses. 
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POWER-OUTPUT  PLOTS- 


SCALE  :  5  . 000E  +  05  UNITS/CM 


Figure  3.6:  Daniel  1  power  spectral  estimates  for  traces  used  in 
the  study. 


see  Cl  owes  (1966,  1969). 

One  final  note  regarding  the  plots  in  figure  3.6, 
is  that  the  power  of  the  filtered  traces  (3  to  55  hertz) 
seems  to  fall  into  two  main  bands,  one  between  10  and  15 
hertz  and  one  between  20  and  30  hertz.  These  bands  likely 
represent  reflected  energy  from  layers  within  the  Earth. 

3.4.  The  Zero-phase  Band-pass  Filter 

The  power  spectra  plots  of  figure  3.6  showed  that 
the  reflected  energy  of  the  seismograms  was  contained 
between  the  frequencies  of  10  to  30  hertz.  To  eliminate 
some  of  the  unwanted  noise  in  the  traces,  for  example 
surface  waves,  the  data  was  subjected  to  a  digital  band¬ 
pass  filter.  During  the  course  of  his  research,  D.  Ganley 
developed  an  eight-pole  zero-phase  shift  recursive  Butter- 
worth  band-pass  filter.  The  Butte rworth  filter  is  design¬ 
ed  to  separate  signal  and  noise  when  they  occur  in  distinct 
frequency  bands.  This  type  of  filter  does  not  exhibit 
Gibb's  phenomena,  a  high  frequency  oscillation  which  super¬ 
imposes  itself  upon  the  amplitude  spectrum  of  the  filter. 
Since  the  filter  is  zero-phase  shift,  there  is  no  lag 
between  different  frequency  components.  Also,  to  prevent 
aliasing  problems,  a  bilinear  Z  transform  was  used  in  the 
design  of  the  filter. 

Basically,  the  program  works  as  follows  (Kanasewich, 
1975).  First,  the  roots  of  the  transfer  function  in  the 
complex  Laplace  transform  variable  (p)  are  calculated. 


From  these  four  poles,  eight  poles  in  complex  conjugate 
pairs,  are  computed  by  a  transformation  of  a  low  to  a  band¬ 
pass  function  (complex  frequency,  s).  These  are  used  to 
calculate  the  filter  gain,  and  the  Z  transform  of  the 
Butterworth  band-pass  filter.  The  latter  is  given  by 

W<Z)  =  B]  (ZJb'uJb'iZJB^Z)  3'4'1 

where  B ^  ,  B  ^  ,  B  ^  and  B  ^  are  the  filter  coefficients 
obtained  from  the  S-poles.  The  filtered  output,  Y ( Z ) ,  is 
obtained  from 


Y(Z)  =  W(Z) • X ( Z)  3.4.2 

where  indicates  a  product.  In  the  discrete  time 

domain,  this  convolution  operator  may  be  expressed  as 

m 

Y+  =  W*X=Y  W.X.  3.4.3 

1  1=0  1  t-1 

where  N  is  the  number  of  points  in  the  input  trace.  A 
recursive  relation  is  used  to  obtain  the  final  output 
(figure  3.7).  To  obtain  zero-phase  shift,  first  the  data 
is  convolved  with  the  filter  to  produce  an  output  trace. 
This  output  trace  is  then  time  reversed  and  again  convolved 
with  the  filter.  The  output  of  this  convolution  is  time 
reversed  to  yield  the  desired  filtered  data.  So  that  the 
gain  of  the  filter  is  unity,  each  point  of  the  output 


X  ( z) 


Input  Data 


*  1  - 


D(z)  - 


where  Di ,  i=l,...,8  are  coefficients  of  z 
and  zaof  the  Butterworth  polynomials. 


Figure  3.7:  A  recursive  band-pass  Butterworth  filter  with 

eight  poles  and  four  zeros  (after  Kanasewich,1975) . 


- 
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trace  is  divided  by  the  gain  of  the  filter. 

In  applying  the  filter  to  the  data,  it  was  decided 
that  a  relatively  wide  pass  band  would  be  used.  Narrow 
band-pass  filtering  makes  the  output  seismogram  quite 
oscillatory,  and  destroys  the  character  of  the  reflected 
pulse.  Therefore,  a  filter  having  cutoff  frequencies  of 
8  and  40  hertz  was  used. 

An  example  of  the  effectiveness  of  the  filter  may  be 
seen  in  figure  3.1.  Here,  we  can  see  that  the  high  fre¬ 
quency  peaks  generated  in  the  clipped  portion  have  been 
totally  eliminated.  However,  this  part  of  the  trace  is 
distorted  and  unusable. 

Figure  3.8  also  exhibits  traces  before  and  after  the 
application  of  the  8  to  40  hertz  filter.  The  first  part 
of  the  figure  shows  a  group  of  traces  that  are  strongly 
disrupted  by  ground  roll.  Beneath  this  group  of  traces 
is  the  filtered  output  showing  no  sign  of  this  long  period 
noise.  Overall,  the  quality  of  the  record  section  can  be 
seen  to  have  been  improved  to  a  large  degree. 

Power  spectral  estimates  were  done  on  the  filtered 
data  (figure  3.9)  to  see  where  the  power  of  the  signal  is 
confined.  As  previously  mentioned,  the  power  is  confined 
to  two  bands,  one  centered  at  13  hertz  and  one  centered  at  27 
hertz.  The  lower  frequency  band  is  probably  the  result 
of  signal  from  deeper  layers,  and  the  higher  frequency 
band  likely  represents  energy  from  shallower  horizons. 


- 


Figure  3.8:  Eleven  traces  of  shot  34B  before  and  after 

8  to  40  hertz  band-pass  filtering.  The 
lower  group  of  traces  are  the  filtered 
output.  The  origin  is  marked  by  the  'O', 
with  the  time  between  adjacent  crosses  being 
1  second. 
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Fi gure  3.9: 


Dan i ell  power  spectral  estimates  for  traces 
after  8  to  40  hertz  band-pass  filtering. 
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3.5. 


Normal  Moveout  Corrections 


Dynamic  corrections  are  displacements  in  time 
between  the  individual  traces  of  a  given  record  that  vary 
with  time.  These  corrections  are  applied  to  simulate  a 
vertical  ray  path  and  so  remove  the  normal  moveout  (or 
time  lag)  which  depends  upon  the  source  to  receiver  dis¬ 
tance  and  the  velocity  distribution  of  the  underlying 
ground  (Mateker,  1965).  From  Fermat's  principle  we  know 
that  the  wavefront  is  propagating  along  the  shortest  time 
paths.  In  regions  where  the  layers  are  horizontal,  or 
have  very  gentle  dips,  the  arrival  time  of  a  reflection 
from  the  nth  layer  at  a  distance  X  is  given  by  the  infinite 
series  (Taner  and  Koehler,  1969) 

T ( x , n ) 2  =  C1  +  C2X2  +  C3X4  +  .. .  3.5.1 


where 


T(x,n)  =  two-way  travel  time  for  the  nth  layer  at 
a  distance  x  from  the  shot  location 
C.j  =  coefficients  dependent  on  interval  velocity 
and  layer  thickness 

A  series  of  model  studies  by  Al-Chalabi  (1973)  have 
demonstrated  the  accuracies  that  are  possible  by  truncat¬ 
ing  equation  3.5.1,  and  have  shown  that  the  accuracy  is 
dependent  upon  the  ratio  of  spread  length  to  reflector 
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depth.  When  this  ratio  was  less  than  1.0,  the  second 
order  approximation  was  correct  to  0.5  percent.  For  a 
ratio  equal  to  2.0,  the  travel  times  were  correct  to  about 
2.0  percent  with  a  two  term  series.  In  our  case,  the  maxi¬ 
mum  shot  to  detector  distance  was  about  10  kilometers, 
yielding  a  ratio  well  below  1.0  for  the  reflectors  of  interest. 
Thus,  a  two  term  truncation  of  the  infinite  series  is  suffi¬ 
cient. 


T(x,n)2  =  C1  +  C2X2 


3.5.2 


Taner  and  Koehler  (1969)  have  calculated  the  coeffic¬ 
ients  C-j  and  C  ^  ,  with  the  result 


T(x  ,n)2  =  T ( o  ,  n ) 2  + 


3.5.3 


where 


T(o,n)  =  two-way  normal  incidence  time  for  the 


nth  1 ayer 


3.5.4 


d(K)  =  thickness  of  the  Kth  layer 


V ( K )  =  interval  velocity  of  the  Kth  layer 


V(n)  =  time  weighted,  root  mean  square  velocity 
as  described  by  Dix  (1955) 


■ 


IK  =  1  V(Krt(K) 
T  (  o  ,  n  ) 


3.5.5 


t ( K )  =  two-way  normal  incidence  time  in  the  Kth 


layer 


3.5.6 


Equation  3.5.3  defines  a  hyperbola  (figure  3.10), 
and  is  accurate  to  within  0.5  percent  of  the  actual  travel 
times  as  given  by  the  exact  solution,  which  is  more  than 
adequate  for  deep  reflections. 

To  remove  the  normal  moveout  in  the  data,  a  two  part 
procedure  was  adopted.  First,  a  model  was  assumed  for  the 
Earth's  crust  in  the  region  of  the  survey.  An  average 
shallow  structure  was  determined  from  a  cross-section  con¬ 
tained  in  Green  et  a  l .  ' s  (1977)  report,  whereas  a  first 
approximation  to  a  deep  structure  is  from  Gurbuz  (1970). 
This  model  is  presented  in  figure  3.11.  Using  this  model, 
a  computer  program  was  written  which  computed  the  travel 
times  for  rays  which  assumed  reflectors  at  intervals  of 
0.04  seconds  of  normal  incidence  time  for  35  shot  to 
detector  distances.  The  travel  time  curves  for  reflectors 
at  the  lower  three  layers  in  the  model  are  shown  in  figure 
3.12,  exhibiting  the  characteristic  hyperbolic  form.  Next, 
a  program  was  written  which  applied  the  normal  moveout 
corrections,  and  it  worked  in  the  following  manner.  First, 
the  program  linearly  interpolated  the  values  for  the  travel 
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Figure  3.10:  Travel  time  curve  for  a  hypothetical  horizontal  reflector. 
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5 

9.59 

7.10 

7.90 

Figure  3.11:  The  crustal  model  used  for  normal  moveout  corrections. 
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Figure  3.12:  Travel  time  curves  for  the  three  deepest  layers 


times  assuming  interfaces  at  0.004  seconds  (the  sampling 
interval)  for  the  35  shot  to  detector  distances.  Next,  for 
each  particular  trace,  the  program  interpolated  an  amplitude 
corresponding  to  the  travel  time  of  a  reflector  at  each 
digital  point.  The  final  result  simulated  a  series  of 
traces  representing  rays  that  appeared  to  have  been  recorded 
with  normal  incidence. 

Since  the  model  did  not  assume  a  continuous  velocity 
increase  with  depth,  some  of  the  early  times  were  double 
valued  for  large  shot  to  receiver  separations.  Therefore, 
in  this  region,  it  was  impossible  to  determine  which  reflec¬ 
tion  contributed  to  the  amplitude  of  the  trace,  and  thus 
this  part  of  the  trace  was  zeroed  (figure  3.13).  In  con¬ 
junction  with  this,  Mair  et  a  Z .  (1976)  point  out  that  off¬ 
sets  should  not  be  so  large  as  to  introduce  large  variations 
in  reflection  character  due  to  varying  reflection  angle. 

For  the  Manitoba  data,  only  a  few  seconds  had  to  be  deleted. 
With  the  Alberta  data  though,  the  high  gain  setting  necessary 
to  record  deep  reflections  resulted  in  almost  five  seconds 
of  clipped  data,  which  I  chose  to  delete  at  this  time. 

Figure  3.14  exhibits  four  traces  before  and  after  the 
removal  of  the  normal  moveout  for  offsets  of  4,  6,  18,  18, 
and  30  shot  to  detector  distances.  The  lower  corrected 
traces  have  been  compressed  throughout,  the  degree  varying 
with  time  along  the  record  and  with  the  shot  to  detector 
spacing.  For  a  spacing  of  four  units,  the  trace  does  not 
change  very  much  after  the  correction  has  been  applied. 


* 
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T(x,n) 


Figure  3.13:  T(x,n)  versus  T(0,n)  for  a  shot  to  receiver  separation 
of  2633  meters  based  on  the  model  of  figure  3.11. 

Times  are  measured  in  seconds. 


' 


Figure  3.14:  Comparison  of  traces  before  and  after  the 

removal  of  normal  moveout.  The  lower  traces 
have  been  corrected.  Note  the  possible 


reflections  beneath  the  'l's. 
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However,  for  a  spacing  of  30  units,  the  early  parts  of  the 
trace  are  greatly  contracted,  with  the  later  sections  being 
only  somewhat  compressed.  At  this  time,  any  primary 
reflections  present  in  the  traces  should  line  up.  Figure 
3.14  suggests  possible  reflections  now  in  almost  perfect 
time  alignment  (beneath  the  '  1  1  s  )  . 


3.6.  Normalization 

With  the  data  being  collected  on  three  different 
systems,  each  with  a  different  gain  setting,  it  was  found 
that  each  trace  had  a  different  amplitude  range.  Thus,  the 
next  step  in  the  analysis  of  the  data  was  the  normalization 
of  the  amplitudes  of  the  traces.  For  this,  it  was  decided 
to  compute  the  average  amplitude,  for  a  particular  time 
window,  of  one  trace  per  group  of  traces  belonging  to  a 
particular  subsurface  point,  and  make  the  rest  of  the  traces 
in  the  group  have  this  average  amplitude.  The  normaliza¬ 
tion  procedure  may  be  expressed  as  follows  (Clowes,  1966): 

n 


Y(i)  =  ( X ( i  )  -  X) 


I1  =  1  |Xs(j)-Xs  | 
I-=1  I  X  ( j  )  -X  | 


3.6.1 


where 


Y  =  normalized  digital  value 
X  =  unnormalized  digital  value 
's'  refers  to  the  standard  trace 

refers  to  an  average  over  the  time  window 


Equation  3.6.1,  as  programmed  by  the  author,  works 
as  follows.  The  average  value  for  the  time  window  was 
calculated  for  all  traces.  This  mean  was  then  removed 
from  all  points  in  the  particular  trace.  Then,  the  sum  of 
the  absolute  value  of  the  deviations  from  the  mean,  within 
the  time  window,  was  computed  for  each  trace.  The  normal¬ 
ized  values  were  obtained  by  multiplying  the  mean  removed 
traces  by  the  ratio  of  the  sum  of  the  deviations  for  the 
standard  trace  to  that  for  each  of  the  other  traces  belong¬ 
ing  to  a  particular  subsurface  point.  This  not  only 
yielded  normalized  traces,  but  within  the  time  window, 
the  D.C.  average  was  removed. 

3.7.  Common  Depth  Point  Stack 

One  of  the  main  problems  encountered  in  reflection 
seismology  is  how  to  remove  the  noise  while  maintaining  a 
high  reflection  signal.  To  overcome  this  problem,  Mayne 
(1962)  introduced  a  multiple  coverage,  common  depth  point 
technique.  This  method  is  a  stacking  process  which  allows 
multiple  coverage  of  the  same  subsurface  points  (Dobrin, 

1976  ;  Telford  et  al .  ,  1976)  utilizing  different  shot  and 
detector  locations  (figure  3.15).  The  accuracy  of  the 
technique  depends  upon  applying  the  correct  time  corrections 
to  each  trace,  both  static  and  dynamic,  so  that  the  primary 
reflections  are  properly  stacked  (Mayne,  1967).  Stacking, 
in  this  manner,  results  in  the  attenuation  of  random  noise, 
multiple  reflections  and  reverberations  (Marr  et  al.  ,  1967; 
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400  percent  coverage. 
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Courtier  et  al . ,  1967). 

A  trace  weighting  technique  for  improving  the  opera¬ 
tion  of  the  common  depth  point  (C.D.P.)  stack  was  suggested 
by  Brown  et  al .  (1977).  This  process  was  designed  for 
vertical  crustal  reflection  data;  data  in  which  the  input 
traces  are  recorded  from  the  same  shot  and  receiver  posi¬ 
tions.  With  different  shot  to  receiver  separations,  as  is 
the  case  for  C.D.P.  data,  the  determination  of  the  weighting 
factors  is  extremely  difficult  due  to  variations  in  reflec¬ 
tion  amplitudes  and  frequencies  between  traces.  White 
(1977)  also  showed  for  low  signal  to  noise  ratios  that  the 
errors  in  determining  the  weighting  factors  were  very  large. 
Therefore,  a  weighted  stack  was  not  used.  Any  traces 
which  were  obviously  contaminated  by  noise  were  simply  not 
included  in  the  stack.  This  would  occur  if  the  seismometer 
output  was  dominated  by  local  traffic,  et  cetera. 

The  average  subsurface  coverage  used  in  the  collec¬ 
tion  of  the  data  was  400  percent,  or  four- fold  (figure  2.6). 
At  this  time,  the  static  and  dynamic  corrections  have  been 
applied  to  the  data,  but  the  static  corrections  may  not 
have  been  exact.  Therefore,  utilizing  the  cross  correla¬ 
tion  program,  all  of  the  traces  for  a  particular  subsurface 
point  were  cross  correlated  together  in  the  vicinity  of 
possible  reflections.  Using  these  coefficients  along  with 
plots  of  the  various  traces,  I  was  able  to  refine  the 
static  time  corrections  for  all  traces.  A  stacking  program 


was  then  written  which  allowed  one  to  shift  the  traces 


according  to  the  revised  statics.  The  operation  of  the 
program  may  be  mathematically  expressed  as 

Y(j)  =  i-  l  X  i  ( j )  3.7.1 

n  i  =1 


where 


X i ( j  )  =  time  corrected  input  traces 
n  =  number  of  input  traces 
Y  (  j  )  =  output  trace 

Figure  3.16  demonstrates  the  effectiveness  of  the 
stacking  process  on  two  types  of  input  signals.  In  both 
cases,  we  observe  an  increase  signal  to  noise  ratio. 

Figures  3.17,  3.18  and  3.19  show  the  output  of  the 
stacking  program  on  three  sections  of  the  data.  In  the 
shallow  sections,  for  subsurface  points  169  to  183,  we  can 
see  coherent  energy  at  about  6.6  and  8.1  seconds.  These 
two  energy  bands  can  be  seen  quite  clearly  on  figure  3.18, 
but  not  very  well  at  all  on  figure  3.19.  In  the  deeper 
sections,  little  bursts  of  coherent  energy  can  be  seen  on 
all  of  the  plots,  notably  around  12.5  seconds.  For  all 
cases,  the  pulses  are  still  quite  obscured  by  the  noise, 
with  the  phase  velocity  of  the  reflections  being  impossibl 
to  determine.  The  bursts  of  coherent  energy  found  at 
various  times  throughout  the  sections  are  most  probably 


Figure  3.16:  Application  of  the  C.D.P.  technique  to  test 


data.  In  (a),  we  have  two  input  channels 
consisting  of  (1,1, 0,0)  and  (1,0, 1,0).  The 
signal  is  contained  in  sample  one.  On  output, 
we  have  (1,0. 5, 0.5,0),  with  a  signal  to 
noise  ratio  of  2,  and  the  noise  being  spread 
out  over  two  samples.  Example  (b)  shows  the 
result  of  a  five- fold  stack  on  a  spike  input 
pulse  and  coherent  noise.  The  noise  arrives 
with  an  apparent  velocity  of  five  points  per 
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Figure  3.17:  Common  depth  stack  for  subsurface  points  169  to 
183.  The  deep  section  (12  to  18  seconds)  is 
amplified  five  times. 
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Figure  3.18:  Common  depth  stack  for  subsurface  points  257  to  267. 
The  deep  section  is  amplified  2.5  times. 
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Figure  3.19:  Common  depth  point  stack  for  subsurface  points  381  to 
391.  The  deep  section  is  amplified  five  times. 


primary  reflections,  based  on  their  amplitudes  and  the  fact 
that  multiples,  being  not  in  proper  time  alignment,  would 
have  been  attenuated  by  the  stack. 

Figure  3.20  exhibits  eleven  traces  corresponding  to 
shot  144,  one  trace  from  each  subsurface  point  from  257  to 
267.  On  the  plot,  coherent  energy  appears  to  exist  at 
about  6.6  seconds  only.  However,  no  distinct  pulse  can  be 
clearly  seen.  In  general,  the  record  is  much  more  noisy, 
and  of  far  less  value  than  the  corresponding  four-fold 
stack  shown  in  figure  3.18. 

To  note  the  frequency  of  the  main  pulse,  we  made  a 
power  spectral  estimate  centered  around  the  pulse  at  6.6 
seconds.  Power  spectral  estimates  were  made  for  each  sub¬ 
surface  point,  and  the  results  were  averaged.  The  time 
window  was  0.2  seconds  in  length,  with  a  taper  being 
applied  to  0.1  seconds  of  the  data.  Several  thousand 
zeros  were  added  to  the  signal  to  increase  the  resolution 
of  the  estimate.  From  figure  3.21,  we  note  that  the  fre¬ 
quency  of  this  main  pulse  is  about  11  hertz,  yielding  a 
wavelength  of  about  590  meters. 

The  common  depth  point  stack  has  been  shown  to  be  a 
powerful  tool  for  enhancing  primary  reflections,  provided 
the  static  and  dynamic  corrections  are  accurate.  Random 
noise  is  attenuated  as  the  square  root  of  the  number  of 
input  channels  (Dobrin,  1976),  while  the  coherent  noise  is 
attenuated  as  a  function  of  its  apparent  velocity.  However, 
as  we  can  see  from  figures  3.18  through  3.20,  primary 
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Figure  3.21:  Average  power  of  the  6.6  second  pulse  for  subsurface  points  257  through  267. 


reflections  are  still  somewhat  obscure,  and  we  would  hope 
to  be  able  to  enhance  the  signals  even  more.  This  will  be 
the  subject  of  the  following  chapter.  In  order  to  apply 
signal  enhancement  techniques  to  the  stacked  data,  the 
stacked  traces  were  normalized  with  respect  to  one 


another 


CHAPTER  4 

MULTICHANNEL  SIGNAL  ENHANCEMENT  TECHNIQUES 

4.1.  The  Velocity  Filter 

Velocity  filters  are  a  class  of  two  dimensional 
linear  operators  which  are  used  to  extract  the  desired 
signals  from  a  background  of  coherent  noise.  The  filter 
operates  in  the  frequency  and  wave  number  space,  and  thus 
requires  a  spatial  array  of  detectors  (Kanasewich,  1975). 

In  the  design  of  the  filter,  it  is  assumed  that  both  the 
signal  and  noise  are  perfectly  coherent  plane-waves,  but 
arrive  at  the  detectors  with  different  apparent  velocities. 
These  filters  have  been  applied  in  exploration  seismology 
for  many  years,  and  have  been  termed  pie-slice,  fan,  velo¬ 
city  or  beam-forming  filters. 

The  operation  of  a  velocity  filter  may  be  stated 
quite  briefly.  By  applying  a  two  dimensional  space-time 
operator  in  the  frequency-spatial  frequency  domain,  all 
events  on  a  multichannel  record,  whose  apparent  velocities 
range  from  -V  to  +V  are  passed.  Thus,  the  desired  transfer 
function  is 


Y(f ,k) 


0 


V 


f  > 


where  K,  the  spatial  frequency,  is  given  by 
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where  A  is  the  wave  number. 

Figure  4.1  illustrates  this  ideal  transfer  function. 

The  impulse  response  of  the  filter  is  determined  by 
the  application  of  the  two  dimensional  inverse  Fourier 
transform  on  the  transfer  function  (Kanasewich,  1975). 
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4.1.3 


The  frequency  limits  are  given  by  the  Nyquist 

limits 


2 AX  <  KN  <  2 AX 

<  f  <  J_ 

2 At  tN  2 At 


4.1.4 


where  AX  is  the  detector  spacing  and  At  is  the  sampling 
interval.  The  time  delay  to  the  nth  sample  point  is 


Tn  =  nAt 


4.1.5 


and,  the  distance  to  the  mth  detector  from  the  center  of 
the  recording  array  is 


Xm 


(m  +  i)AX 


4.1.6 
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k=spatial  frequency 
kM=Nyquist  spatial  frequency 

f=frequency 
fN=Nyquist  frequency 


In  the  pass  region,  Y(f,k)=l 
In  the  reject  region,  Y(f,k)=0 


Figure  4.1:  The  velocity  filter  transfer  function. 


Utilizing  equations  4.1.3  through  4.1.6,  Embree  et  al . 
(1963)  have  shown  that  the  impulse  response  of  the  filter 
may  be  expressed  as 


W(Tn  ,Xm) 


4.1.7 


This  impulse  response  is  time  and  space  symmetric, 
has  zero-phase  shift,  and  may  be  convolved  with  the  origi¬ 
nal  traces  to  produce  a  velocity  filtered  record  (Kanase- 
wich,  1975).  However,  this  classical  method  is  time 
consuming,  as  for  N  input  traces,  N/2  convolutions  must  be 
performed  to  produce  an  output  trace.  To  reduce  the 
computational  time,  Treitel  et  al .  (1967)  transformed  the 
filter  equations  into  the  Z  domain,  and  then  made  use  of  a 
recursive  relation  (Shanks,  1967)  so  that  only  one  convolu¬ 
tion  is  required  for  each  output  trace.  Figure  4.2  depicts 
this  algorithm,  and  is  explained  as  follows.  Each  input 
trace  from  a  given  record  containing  an  even  number  of 
equally  spaced  traces  is  subjected  to  two  time  shifts. 

From  this  pair  of  traces,  a  difference  trace  is  evaluated, 
and  the  appropriately  weighted  difference  traces  are  then 
summed.  The  resulting  single  output  trace  is  then  convolved 
with  the  filter  operator  producing  the  desired  output  trace. 
This  algorithm  is  easily  programmed,  and  computationally 
far  less  time  consuming. 

During  the  course  of  his  Ph.D.  thesis  at  the 
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Figure  4.2:  A  four  channel  velocity  filter  (after  Kanasewich ,1975) . 


University  of  Alberta,  Clowes  (1969)  developed  a  velocity 
filter  program  (appendix  3)  which  allowed  a  filter  with  an 
even  number  of  traces.  Due  to  the  subsurface  coverage 
being  anywhere  from  300  to  700  percent,  it  was  decided 
that  the  application  of  a  four  trace  fan-pass  filter  would 
be  sufficient. 

To  demonstrate  the  effectiveness  of  the  fan-pass 
filter,  seven  different  dipping  events  with  apparent  velo¬ 
cities  ranging  from  12.2  kilometers  persecond  through  an 
infinite  apparent  velocity  to  -12.2  kilometers  per  second 
were  constructed  in  the  input  channels.  These  are  shown  in 
the  upper  part  of  figure  4.3.  The  three  sets  of  traces 
below  the  input  channels  are  the  result  of  applying  three 
different  pass  bands  to  the  input  data.  From  the  diagram, 
it  is  obvious  that  events  with  apparent  velocities  within 
the  specified  range  are  passed  without  distortion,  while 
those  outside  of  the  range  are  greatly  attenuated.  An 
event  which  propagates  at  one  of  the  cutoff  velocities  is 
attenuated  in  amplitude  by  approximately  one-half. 

The  velocity  filter  was  then  applied  to  the  common 
depth  point  stacked  data  of  figures  3.17  to  3.20.  Velocity 
filtered  seismograms  for  these  data  sets  are  contained  in 
figures  4.4  to  4.7.  In  the  shallow  section  of  figure  4.4, 
a  very  strong  reflection  can  be  seen  at  about  6.6  seconds 
(A).  This  event  can  also  be  seen  on  figures  4.5  (B)  and 
4.6  (A).  From  the  plots,  this  reflection  appears  to  result 
from  a  horizontal  layer  in  the  Earth  in  figure  4.4  and  4.5, 
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Figure  4.3:  Application  of  an  eight  channel  velocity 


filter  to  synthetic  data.  Three  different 
pass  bands  are  applied  to  the  data. 
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Figure  4.4:  Application  of  a  four  channel  velocity 

filter  to  the  stacked  data  for  subsurface 
points  169  to  183.  The  deeper  section  (12  to 
18  seconds)  is  amplified  five  times.  Coherent 
energy  may  be  seen  above  A,  B,  C,  and  D. 
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Figure  4.5:  Application  of  a  four  channel  velocity  filter  to  the 
stacked  data  for  subsurface  points  257  to  267.  The 
deep  section  is  amplified  2.5  times. 
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Figure  4.6:  A  four  channel  velocity  filter  as  applied  to  the 
stacked  data  for  subsurface  points  381  to  391.  The 
deep  section  is  amplified  five  times. 


and  from  a  westerly  dipping  layer  in  figure  4.6.  A 
westerly  dipping  event  at  about  8.1  seconds  is  strongly 
seen  on  figures  4.4  (B)  and  4.5  (A),  but  is  somewhat 
obscure  on  figure  4.6  (B).  Other  coherent  phases  of  lesser 
amplitude  are  seen  throughout  the  filtered  sections,  but 
these  are  likely  due  to  random  noise  producing  a  coherent 
line  up  of  energy,  or  is  what  is  left  of  multiple  reflec¬ 
tions  after  attenuation  in  the  C.D.P.  stack. 

In  the  deeper  regions  of  the  velocity  filtered 
sections,  coherent  energy  can  be  seen  at  about  12.5  seconds 
(pulse  C)  on  all  of  figures  4.4  through  4.6.  From  the  pie- 
slice,  this  event  is  dipping  to  the  east  on  figures  4.5 
and  4.6,  and  to  the  west  on  figure  4.4.  Other  deeper  events 
appear  to  exist  on  the  section,  notably  around  14  seconds 
(  pul se  D)  . 

Figure  4.7  exhibits  the  results  of  velocity  filter¬ 
ing  the  traces  belonging  to  shot  144  (figure  3.20),  that  is, 
using  only  one-fold  coverage.  The  traces  were  taken  from 
subsurface  points  257  to  267,  whose  velocity  filtered 
section  may  be  seen  in  the  first  half  of  figure  4.5.  The 
good  quality  reflections  at  6.6  (B)  seconds  is  still 
enhanced,  while  the  dipping  event  at  8.1  (A)  seconds  is 
only  partially  brought  out.  Other  coherent  phases  can  be 
seen  at  9.3  and  11.6  seconds.  These  events  can  also  be 
located  on  figure  4.5,  but  they  are  of  far  less  amplitude. 
Thus,  it  is  believed  that  these  two  events  are  the  result 
of  multiply  reflected  energy  from  a  shallower  horizon. 
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Figure  4.7:  The  velocity  filtered  section  for  the  data  of  figure 

3.20.  Reflected  events  are  seen  at  A  and  B,  with  possible 
multiple  reflections  beneath  the  two  1  I ' s . 
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When  stacked  with  other  traces  of  varying  shot  to  detector 
distances,  they  are  not  quite  in  time  alignment,  and  are 
thus  severely  attenuated.  However,  when  taken  alone, 
multiple  reflections  along  adjacent  traces  would  appear  to 
have  arrived  with  some  apparent  velocity  across  the  spread. 
Therefore,  when  we  velocity  filtered  the  data  section,  the 
multiple  reflections  were  enhanced,  and  are  of  relatively 
high  amplitude. 

The  fan-pass  filter  can  clearly  be  seen  to  have 
improved  the  quality  of  the  reflection  seismograms.  Reflec¬ 
ted  events  occurring  around  6.6,  8.1  and  12.5  seconds, 
which  had  been  obscured  by  noise  on  the  C.D.P.  stacks,  were 
enhanced  greatly  on  all  three  sections.  With  this  techni¬ 
que,  not  only  are  the  primary  reflections  enhanced,  but 
their  apparent  velocities,  and  thus  the  dips  of  the  horizons 
are  easily  determined.  Therefore,  through  the  use  of  this 
technique,  it  should  be  fairly  easy  to  detail  the  crustal 
structure  of  a  particular  region. 

4.2  The  N-th  Root  stack 

A  method  of  non-linear  velocity  filtering  was 
developed  by  Muirhead  (1968),  and  expanded  upon  by 
Kanasewich  et  al .  (1973),  which  they  termed  the  N-th  root 
stack.  Basically,  the  filter  involves  the  extraction  of 
the  N-th  root  of  each  element  in  all  of  the  traces,  where 
N  is  any  positive  integer.  The  traces  are  then  summed  and 
raised  back  up  to  the  N-th  power.  The  filter  was  shown  to 


be  effective  in  removing  random  noise,  with  the  effective¬ 
ness  increasing  with  the  number  of  channels.  For  a  delta 

function  noise  spike,  the  rejection  varies  approximately 
_  n 

as  K  ,  where  K  is  the  number  of  channels.  Unlike  linear 
filters,  one  cannot  obtain  an  impulse  response  or  a 
transfer  function.  Also,  the  output  of  the  filter  is  dis¬ 
torted  to  some  degree. 

The  N-th  root  stacking  filter  may  be  mathematically 
expressed  as 


Yi  v 


Ri  v  Ri  v 


4.2.1 


where 


Ri  v 


K 


1 

IK  w(j) 
j=l 


M'+T.  )  ,  j 
I x( i +t j ) ,  j 


GW(,j) 

“GUT 


I  X  ( i  +T  j  )  ,  j  I 
4.2.2 


with 


Yiv  =  single  channel  output 

X  ( i  )  ,  j  =  ith  data  point  on  the  jth  channel 

N  =  any  positive  integer 

K  =  number  of  channels 

G  =  common  gain  for  normalization 

G  ( j )  =  gain  of  the  jth  channel 

W  (j)  =  weighting  coefficient  for  the  jth  channel 


Also,  we  have 
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where 


D.  =  distance  from  array  center  to  jth  detector 

J 

V  =  apparent  velocity  of  the  wave 

In  the  case  of  normalized  common  depth  point  data, 
these  equations  reduce  to 
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4.2.4 

Thus,  we  are  attempting  to  enhance  particular 
reflections  with  particular  apparent  velocities.  If  N  =  1  , 
we  obtain  a  linear  'delay  and  sum'  velocity  filter. 

Figure  4.8  demonstrates  the  operation  of  the  filter. 
The  upper  part  of  the  figure  consists  of  three  examples, 
each  with  two  input  channels.  One  of  the  channels 
(channel  1)  contains  only  signal,  whereas  channel  2  contains 
both  signal  and  noise.  For  only  pure  signal  (case  a),  it 
is  passed  without  distortion.  In  case  b,  where  only  noise 
is  present,  the  effect  of  increasing  the  power  of  the  stack 
is  evident.  Example  c  has  90  units  of  noise  added  to  the 
signal  to  make  up  channel  2,  and  again  we  note  the  effect 
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Figure  4.8:  The  N-th  stack  as  applied  to  spike  model  data. 


of  increasing  the  power  of  the  stack  in  reducing  the  noise. 

Figure  4.9  demonstrates  the  effect  of  the  filter  on 
a  synthetic  pulse.  The  input  consists  of  five  channels 
each  comprising  two  input  signals.  The  signals  are  some¬ 
what  distorted  due  to  the  effect  of  a  slanting  noise  field 
arriving  with  an  apparent  velocity  of  five  digital  points 
per  channel.  As  we  can  see,  the  two  input  signals  are 
reproduced,  with  some  distortion,  with  the  noise  field 
being  totally  eliminated.  This  indicates  that  this  filter 
may  be  used  to  determine  the  apparent  velocities  of  the 
arrivals,  and  thus  the  dips  of  the  underlying  layers. 

A  program  was  written  which  would  perform  an  eight 
root  stack  (appendix  3)  using  the  normalized  common  depth 
point  data  as  input  (shallow  sections).  The  number  of 
channels  to  use  in  the  stack  was  chosen  to  be  four  (the 
same  as  used  in  the  fan-filter),  and  stepouts  were  chosen 
to  range  from  -4  to  +4  digital  points.  The  output  of  this 
program  on  the  C.D.P.  traces  169  to  183  is  seen  in  figure 
4.10.  We  see  that  strong  arrivals  occur  at  stepouts  equal 
to  zero  and  three,  for  reflections  at  6.6  (A)  and  8.1  (B) 

seconds.  The  signals  are  a  bit  distorted,  but  the  noise 
is  so  well  destroyed  that  the  exact  timing  of  the  reflec¬ 
tion  is  possible.  For  stepout  equal  to  zero,  the  slight 
jump  of  the  pulse  between  output  channels  8  and  9  is  very 
noticeable,  and  very  accurately  timed.  This  jump  can  also 
be  seen  on  the  velocity  filter  plots  (figure  3.17),  but 
not  quite  so  clearly. 
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Figure  4.9:  Application  of  an  eighth  root  stack  to  a  spike  input  pul 
and  coherent  noise.  The  noise  arrives  with  an 
apparent  velocity  of  five  points  per  trace. 


. 


' 


Figure  4.10:  An  eight  root  stack  of  the  stacked  data  for 


subsurface  points  169  to  183. 
energy  may  be  seen  above  A  and 
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The  output  of  the  filter  for  subsurface  points  257 
to  267  may  be  seen  in  figure  4.11,  and  may  be  compared  with 
the  velocity  filtered  output  seen  in  figure  3.18.  Again, 
for  zero  stepout  we  have  a  strong  arrival  at  6.6  seconds 
(A);  and  for  a  stepout  of  three,  we  have  a  strong  arrival 
at  8.1  seconds  (B).  These  events  are  quite  clear  in 
figure  3.18,  but  the  remaining  noise  which  is  also  seen  in 
figure  3.18  is  much  more  reduced  by  the  eight  root  stack. 

For  subsurface  points  381  to  391,  the  action  of  the 
filter  may  be  seen  in  figure  4.12.  With  stepout  equal  to 
three,  we  have  coherent  energy  at  6.7  (A)  and  8.3  (B) 
seconds.  The  pulse  at  6.7  seconds  is  quite  clear  on  figure 
3.19,  but  the  event  at  8.3  seconds  is  barely  noticeable. 

The  N-th  root  stack,  multichannel  filter  is  very 
useful  in  enhancing  seismic  reflection  data.  The  effective 
ness  of  the  filter  increases  with  increasing  power  of  the 
stack,  and  the  filter  effectively  removes  both  random  noise 
and  noise  arriving  with  finite  apparent  velocity.  In  the 
case  of  the  fan-pass  filter,  the  apparent  velocity  of  the 
arriving  wave  may  be  determined  within  a  certain  pass  band. 
However,  with  the  N-th  root  stacking  filter,  the  apparent 
velocity  of  the  arrivals  may  be  determined  exactly.  This 
great  resolution  is  very  important  when  dealing  with  reflec 
tions  from  dipping  beds. 

4.3.  Multichannel  Adaptive  Filter  Theory 

The  two  main  types  of  adaptive  multichannel  filters 
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Figure  4.11:  The  eight  root  stack  as  applied  to  the  stacked  data  for 
subsurface  points  257  to  267. 
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Figure  4.12:  Application  of  an  eight  root  stacking  filter  on  the 
stacked  data  for  subsurface  points  381  to  391. 


applied  to  reflection  seismic  information  are  the  maximum 
likelihood  filter  and  the  Wiener  filter.  Adaptive  multi¬ 
channel  filters  are  constructed  from  the  input  data,  which 
is  often  common  depth  point  data.  Utilizing  the  autocor¬ 
relation  and  cross  correlation  statistics  between  the 
input  traces,  the  Wiener  filter  compresses  the  unknown 
signal  to  a  delta  function  (Mercado,  1973).  Therefore, 
primary  reflected  energy  will  be  noted  by  a  significant 
peak  in  the  zero  lag  cross  correlation  term.  On  the  other 
hand,  maximum  likelihood  filters  are  constructed  to  detect 
coherent  signal  arriving  with  infinite  apparent  velocity 
across  the  spread,  and  to  minimize  the  output  power  of  the 
traces.  Random  noise  and  noise  arriving  with  finite 
apparent  velocity  are  thus  reduced. 

Maximum  likelihood  filters  were  developed  to  analyze 
data  obtained  from  nuclear  blast  detectors  in  Arizona  and 
Montana  (Mercado,  1978).  The  multichannel  data  is  first 
studied  for  event  free  time  windows  near  intervals  of  sig¬ 
nal.  Through  various  digital  techniques,  the  coherent 
signals  are  adjusted  to  have  the  proper  time  alignment. 
Then,  using  the  noise  interval,  the  filter  is  designed, 
and  subsequently  applied  to  a  time  window  containing  both 
the  signal  and  the  noise. 

Multichannel  filtering  is  a  process  that  makes  use 
of  data  recorded  digitally  on  K  channels,  X^(t),...,X^(t), 
where  each  X .  ( t )  represents  one  channel  of  seismic  data. 

A  digital  seismic  trace  may  be  represented  by  a  finite 


' 


series  of  real  numbers  which  are  samples  of  the  seismic 
energy  being  recorded  at  a  constant  time  interval  At. 


1  1 1 


The  autocorrelation  of  each  X.(t)  is  defined  to  be 

rT  =  ^  { X i ( t ) • X . (t-T )T}  t  =  0,±1,±...  4.3.1 

where  N  is  the  number  of  samples  of  X . ( t )  and  indicates 

inner  product  and  '  T '  means  transpose. 

The  cross  correlation  between  any  two  seismic 

channels  X.(t)  and  X.(t)  is  defined  to  be 
i  J 

gT  =  1  {X . ( t) -Xj ( t-x)T}  t=0,±1,±...  4.3.2 

In  the  multichannel  case,  we  can  represent  the  multichannel 
data  by  a  column  vector  X(t)  (Wiggins  et  al . }  1965) 

X ( t )  =  (X1  ( t) ,  .  . . ,XK(t) )T.  4.3.3 

The  multichannel  autocorrelation  function  now 

becomes 


rT  =  ^  {X(t)*X(t-x)T}  t  =  0  ,  ±1  ,  ±  .  .  .  4.3.4 

where  '*'  indicates  matrix  multiplication.  These  coeffi¬ 
cients  satis fy  (Robinson,  1967) 

T 


r-x 


4.3.5 
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The  multichannel  filter  is  represented  by  the 


coefficients  f ( 0 f ( t ) ,  where  each  f(i)  is  a  Kxl 
matrix.  The  multichannel  signal  is  the  input  to  the  filter 
with  the  resulting  single  channel  output,  0(t).  To  obtain 
the  output,  one  must  convolve  the  input,  X(t),  with  the 
filters. 


0(t)  =  X ( t ) I  f(0)  +  ...  +  X(t-T)J  f ( t )  . 


4.3.6 


Using  appendix  1,  equation  4.3.6  may  be  rewritten 


a  s 


0  (  t  )  =  X.J  (  t-T  )  *fi  (  T  )  +  ...  +  X  K  (  t-T  )  *f  K  (  T  )  4.3.7 


where  each  f  .  ( t )  is  a  (  t  + 1  )  x  1  matrix. 

If  the  filter  is  to  pass  coherent  signal  between 


input  channels  undistorted,  then  the  sum  of  the  filters 
must  be  a  delta  function  (Mercado,  1978). 


f  1  (  t  )  +  ...  +  f  k(t  )  =  <5  (x  ) 


4.3.8 


or 


f-|  (  t  )  =  6  (  t  )  -  f  £  (  r  ) 


4.3.9 


Substituting  4.3.9  into  4.3.7, 


1  1  3 


0(t)  =  X.|  (  t-T  )  *{  6  (  T  )  -  f  2  (  T  )  -  .  .  .  -f  K  (t  )  }  +  .  .  .  +XK  (  t-T  )  *f  k(t  ) 


=  X1  (t)-X1  (t-T)*f2(T)- - -X1  (t-T)*f k(t)  +  .  .  . 


+XK(t-x)*fK(T) 

=  X]  (t)+[X2(t-x)-X1  (t-x)]*f2(x)+..  . 

+  [XK(t-x)-X1 ( t-x ) ] *f K ( x )  .  4.3.10 

Therefore , 

0 ( t )  =  X ^ ( t ) +A^ ( t -x ) *f 2 ( x ) + . . . +A^ ( t -x ) *f ^ ( x )  4.3.11 

where  A . ( t - x )  =  X  .  ( t - x )  -  X^(t-x). 

The  output  power  is  given  by  0(t)"l"*0(t),  and  the 
minimization  of  this  power  will  yield  the  least-squares 
solution  for  f 2 ( x ),..., f^ ( x )  (see  appendix  2  for  details 
regarding  the  three  channel  case).  The  resulting  simul¬ 
taneous  equations  necessary  to  obtain  the  filter  coeffi¬ 
cients  are  of  the  form, 
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where  the  first  matrix  is  the  multichannel  autocorrelation 
matrix  of  the  differences  traces  A2(t-x),...,A^(t-x). 

From  equation  4.3.5,  we  note  that  this  matrix  is  Toeplitz 
in  form.  Each  of  the  elements,  r  ,  is  a  (K-l)  x  (K-l) 
matrix.  The  matrix  elements  f(o),...,f(x)  are  (K-l)  x  1 
matrices  representing  the  filter,  and  the  matrix  elements 
g  ,...,g  are  (K-l)  x  1  cross  correlation  terms  between  the 
first  channel,  X ^ ( t ) ,  and  (A2(t-x),...,AK(t-x)}. 

The  solution  of  this  set  of  simultaneous  equations 
gives  the  coefficients  f(o),...,f(x),  which  are  related  to 
f  2  (  t  )  ,  ••  •  »  f  K  (  t  )  by 

f( 3 )  =  f2(j ) » • • • >fK^ j=0,...,x  4.3.13 

Thus,  using  equation  4.3.11,  we  are  able  to  obtain 
a  filtered  output  channel.  A  three  channel  maximum  likeli¬ 
hood  filter  program  was  written  based  on  appendix  2. 

Figure  4.13  demonstrates  the  action  of  a  two  channel 
maximum  likelihood  filter  on  a  simple  spike  model.  A 
simple  stack  of  the  channels  yields  a  2  to  1  signal  to 
noise  ratio  (SNR).  If  we  apply  a  three  filter  operator  to 
the  channels,  as  depicted  in  the  figure,  we  obtain  an 
output  channel  with  a  4  to  1  signal  to  noise  ratio.  We 
also  note  that  the  noise  has  now  been  spread  over  four 
sampl es  . 

At  first,  an  attempt  was  made  to  calculate  f(o), 

.  .  .  ,  f  ( x  )  based  on  an  algorithm  of  Wiggins  (1965)  which  makes 
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Figure  4.13:  The  two  channel  maximum  likelihood  filter  as 
applied  to  a  spike  model. 
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use  of  the  multichannel  autocorrelation  matrix  being  in 
Toeplitz  form.  Utilizing  an  inversion  scheme,  based  on 
this  algorithm,  found  in  Robinson  (1967),  a  program  was 
written  to  do  maximum  likelihood  filtering.  However,  this 
technique  did  not  work  at  all,  and  is  likely  the  result  of 
many  of  the  columns  in  the  autocorrelation  matrix  nearly 
being  equivalent  (degenerate  matrix).  This  fact  was  also 
noted  by  Mercado  (1978).  Attempts  at  writing  higher  degree 
maximum  likelihood  filter  programs  were  also  made  (four 
and  five  channels),  but  these  also  met  with  failure.  Again, 
this  is  likely  a  result  of  near  dependency  amongst  the 
columns  of  the  autocorrelation  matrix. 

Thus,  it  was  found  that  a  three  channel  maximum 
likelihood  filter  was  best  programmed.  The  algorithm  for 
solving  the  simultaneous  equations  of  equation  4.3.12  is 
contained  in  the  IMSL  library  of  the  University  of 
Alberta's  computing  center,  and  the  complete  program  is 
located  in  appendix  3. 

Figure  4.14  shows  a  test  function  for  the  program, 
and  the  resulting  output.  As  we  can  see,  the  slanting  noise 
events  are  almost  totally  annihilated,  whereas  the  coherent 
signal  is  well  enhanced  and  accurately  reproduced.  The 
filter  operator  is  40  points  in  length. 

The  application  of  the  40  points,  three  channel 
maximum  likelihood  filter  on  real  data  may  be  seen  in 
figure  4.15.  For  this  section,  I  chose  the  eleven  sub¬ 
surface  points  169  to  179,  whose  C.D.P.  stack  may  be  seen 
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3— CHANNEL  FILTER. 

DISTANCE  BETWEEN  TRACES  IS  0.146  KM 

TIME  IN  SECONDS 

6.40  6.69  6.78  6.97  7.16  7.34 

! - 1 - 1_ I - 1 - 1 

INPUT  PULSE 

Ar 
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Figure  4.14:  A  three  channel  maximum  likelihood  filter  applied  to 

a  spike  input  pulse  and  coherent  noise.  The  noise  arrives 
with  an  apparent  velocity  of  three  points  per  trace.  The 
output  is  amplified  five  times. 
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I  — .  i  II  .  J 

6.40  Time  (sec)  7. 40 

Figure  4.15:  Application  of  a  three  channel  maximum  likelihood  filter 
on  the  data  for  subsurface  points  169  to  183.  Coherent 
energy  may  be  seen  beneath  the  'I'  in  the  upper  half. 
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in  figure  3.17  (shallow  section).  The  top  part  of  the 
figure  exhibits  the  output  after  one  pass  with  the  filter, 
and  the  lower  part  exhibits  the  output  after  another  pass 
with  the  filter.  In  the  figure,  we  see  that  the  noise  is 
quite  well  reduced.  The  coherent  pulse  at  6.6  seconds  has 
been  enhanced  somewhat,  but  it  also  has  been  distorted. 

The  slight  step  in  the  pulse  is  also  noticeable.  A  third 
pass  of  the  filter  on  the  data  resulted  in  a  depreciation 
of  the  output . 

The  filter  was  also  applied  to  the  stacked  data 
corresponding  to  the  shallow  section  for  subsurface  points 
257  to  267  (figure  3.18),  with  the  result  being  displayed 
in  figure  4.16.  A  highly  distorted  pulse  may  be  seen 
around  6.6  seconds  which  corresponds  to  the  input  pulse. 

On  the  whole,  the  output  is  not  very  good,  and  further 
application  of  the  filter  to  the  data  did  not  improve  it 
any  more. 

Application  of  the  maximum  likelihood  filter  to  real 
seismic  data  was  not  deemed  very  successful.  Reflecting 
events  were  not  enhanced  very  much,  and  in  some  cases, 
were  distorted  beyond  recognition. 
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Figure  4.16:  Application  of  a  three  channel  maximum  likelihood  filte 
on  the  stacked  data  for  subsurface  points  257  to  267. 
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CHAPTER  5 
CONCLUSIONS 


Seismograms  recorded  in  southeastern  Saskatchewan 
and  southwestern  Manitoba  clearly  indicate  the  existence 
of  near  vertical  incidence,  deep  reflecting  events  from 
within  the  Earth's  crust.  Our  own  reflection  profile 
was  obtained  along  three  lines  in  southwestern  Manitoba 
over  a  distance  of  about  7  2  kilometers.  The  data  was 
recorded  digitally  on  magnetic  tape.  Through  the  use  of 
the  University  of  Alberta's  Amdahl  470V/6  computer,  static 
and  dynamic  corrections  were  applied  to  the  data.  Power 
spectral  analysis  of  the  various  traces  indicates  that  the 
reflected  energy  is  contained  between  the  frequencies  of 
10  and  30  hertz.  Therefore,  a  zero-phase  band-pass  filter 
with  cutoff  frequencies  of  8  and  40  hertz  was  applied  to 
the  data  to  improve  its  overall  quality.  A  major  reflection 
recorded  at  a  two-way  travel  time  of  about  6.6  seconds  was 
found  to  have  a  dominant  frequency  of  about  11  hertz. 

Using  cross  correlation  curves  and  various  plots, 
primary  reflections  were  brought  into  proper  time  alignment 
for  three  sections  of  subsurface  points.  Traces  for  each 
subsurface  point  were  subsequently  stacked  and  investigated 
for  coherent  energy.  To  further  enhance  any  primary 
reflections,  the  techniques  of  velocity  filtering,  N-th 
root  stacking  and  maximum  likelihood  filtering  were  applied 
to  the  data.  Both  velocity  filtering  and  N-th  root  stack¬ 
ing  enhanced  well  all  reflections,  with  maximum  likelihood 
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filtering  having  very  little  success.  To  produce  eight 
velocity  filtered  traces  (each  six  seconds  in  length)  from 
eleven  input  channels  on  the  Amdahl  470V/6  took  approxi¬ 
mately  4.8  seconds,  while  the  corresponding  eight  root 
stack  time  was  about  1.5  seconds. 

Maximum  likelihood  filtering  was  not  found  to  be  a 
very  useful  technique.  Accurate  inversion  of  the  large 
autocorrelation  matrices  is  very  difficult  and  very 
expensive.  To  compute  nine  output  traces  (each  one  second 
in  length)  from  eleven  input  channels  requires  roughly 
13.5  seconds  of  computer  time.  Also,  it  was  found  that 
the  filter  is  very  unstable  when  the  noise  is  not  coherent. 
Thus,  this  technique  should  not  be  used  in  the  analysis  of 
data  where  random  noise  predominates. 

This  thesis  illustrates  that  the  reflection  method, 
with  its  increased  resolution,  can  add  substantially  to 
the  knowledge  of  the  crust.  Suggestions  for  future  studies 
of  this  type  include: 

(1)  Use  of  tapered  arrays  and  multiple  holes  to 
to  increase  the  quality  of  the  reflection 
seismograms.  Also,  a  better  site  selection 

for  the  geophones  and  stringers  (not  near  power 
lines  or  on  gravel  roads)  would  increase  the 
quality  of  the  records. 

(2)  Multiple  reflections  did  not  seem  to  be  a 
problem,  but  were  noticed  on  some  sections.  The 
increase  from  400  to  1200  percent  coverage 
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would  probably  eliminate  these  and  enhance  the 
continuity  of  primary  reflected  energy. 

(3)  A  velocity  analysis  could  be  used  to  determine 
exactly  the  velocities  necessary  to  remove  the 
normal  moveout  (NMO)  associated  with  the  records. 
This  is  not  a  vital  step  in  the  analysis  since 
even  for  spread  lengths  of  up  to  ten  kilometers, 
the  effect  of  variations  in  velocities  on  NMO, 
for  long  ray  paths,  is  very  small.  However,  it 
would  simplify  the  steps  necessary  to  bring 
primary  reflections  into  proper  time  alignment. 

(4)  Band-pass  and  velocity  filtering  can  be  used 
before  stacking  to  remove  various  forms  of 
noise  (for  example,  ground  roll).  With  a 
twelve- fold  stack,  the  multiple  reflections 
should  be  well  attenuated  eliminating  any  need 
for  deconvolution.  An  eight  root  stack  subse¬ 
quently  applied  to  the  stacked  data  should  now 
well  enhance  any  primary  reflections. 

On  the  velocity  filtered  sections  obtained  in  this 
study,  three  main  reflecting  horizons  can  be  seen;  one  at 
about  6.6  seconds,  one  at  8.1  seconds,  and  one  at  12.5 
seconds.  An  approximate  average  velocity  down  to  these 
horizons  may  be  determined  by  using  the  structure  model  of 
figure  3.12,  and  the  following  equation  from  Taner  and 
Koehler  (1969). 
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where  the  parameters  are  the  same  as  those  of  secton  3.5. 
Using  equation  5.1.1,  we  get  average  velocities  to  the 
three  layers  of  about  5.7,  5.9  and  6.2  kilometers  per 
second  respectively.  From  figure  5.1,  the  approximate 
dip  angle  to  each  of  these  horizons  is 
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Figures  5.2  through  5.4  exhibit  the  resultant 
migrated  crustal  sections  with  the  approximate  dips  shown. 
In  order  to  confirm  the  continuity  of  these  horizons,  they 
must  be  located  on  many  more  sections  along  the  three 
lines.  If  the  reflections  are  from  a  transition  zone 
instead  of  a  sharp  velocity  discontinuity,  the  events  may 
not  correlate  over  great  distances  (see  section  1.2). 

These  models  can  be  compared  with  the  preliminary 
refraction  model  for  the  east  side  of  the  Churchill- 
Superior  boundary  (figure  5.5).  The  shallow  layers  seen 
in  the  refraction  profile  were  not  observed  due  to  limita¬ 
tions  in  our  recording  instrumentation.  The  upper  two 
layers  in  the  reflection  models  were  not  detected  by  the 
refraction  method,  possibly  due  to  wide  spacing  in  the  shot 
to  receiver  locations,  but  our  deep  reflector  at  12.5 
seconds  coincides  quite  closely  with  the  third  layer,  and 
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A  B 


dip  of  the  reflector 

Tl=  two  way  travel  time  to  the  reflector  below  A 
T2=  two  way  travel  time  to  the  reflector  below  B 
Therefore, 

Yl=  distance  travelled  by  the  wave  below  A 

=Vav  *  T1 
2 

Y2=  distance  travelled  by  the  wave  below  B 

=Vav  *  T2 
2 

Y3=  Y2-Y1 

=_Vav_*  (T2-T1 )  =  Vav  *AT 
2  2 

Figure  5.1:  The  expression  for  the  dip  angle  for  a  reflector 
where  the  average  velocity  down  to  the  layer  is 
constant. 
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Figure  5.2:  Deep  crustal  model  as  derived  from  subsurface  points  169  to  183. 
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Figure  5.3:  Deep  crustal  model  as  derived  from  subsurface  points  257  to  267. 
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Figure  5.4:  Deep  crustal  model  as  derived  from  subsurface  points  381  to  391. 
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Figure  5.5:  Preliminary  findings  of  the  east-west  refraction  profile  (after  Green, 1977). 


the  energy  at  about  14.0  seconds  coincides  approximately 
with  the  Mohorovicic  discontinuity  seen  in  the  refraction 
profile.  The  refraction  profile  across  the  boundary  zone 
seems  to  indicate  a  boundary  at  100.9  degrees  longitude 
corresponding  to  the  break  in  the  second  layer.  Across 
this  boundary,  the  Mohorovicic  discontinuity  is  at  a  depth 
of  about  46  kilometers.  Unfortunately,  this  structure 
break  occurs  at  the  very  end  of  our  reflection  profile, 
and  thus  in  order  to  see  it,  the  reflection  profile  should 
be  repeated  over  this  region.  Since  the  field  program 
appears  likely  to  be  funded  in  1979,  it  is  recommended  that 
this  interesting  area  be  studied  in  detail. 

In  conclusion,  the  velocity  filtered  sections  clearly 
enhanced  all  reflections  within  a  specified  pass  band,  but 
only  limits  could  be  placed  on  the  phase  velocity.  With 
the  N-th  root  stack,  only  events  with  the  proper  time 
alignment  for  a  specified  stepout  were  passed,  with  all 
other  events  being  eliminated.  The  technique  is  very 
dependent  upon  amplitude  and  proper  time  alignment  of  the 
various  traces,  which  when  all  are  correct,  make  accurate 
computations  of  phase  velocity  very  simple.  The  N-th  root 
stack  has  been  successfully  used  to  enhance  low  frequency 
refraction  seismic  records  (Kanasewich  et  al .,  1973). 
However,  this  thesis  is  the  first  demonstration  of  the 
value  of  the  N-th  root  stack  in  the  signal  enhancement  of 
reflection  seismic  profiles. 
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APPENDIX  1 

MULTICHANNEL  FILTER  OUTPUT  EQUATIONS 


The  output  of  a  multichannel  filter  may  be  expressed 
as  (equation  4.3.6) 

0 ( t )  =  X(t)T*f(o)+. . .+X(t-T)T*f (t) 


X-|  ( o )  ...  XK(o) 

f-l  (o) 

X-j  (-x)  . . .  XK(-x) 

f-|  (t) 

•  • 

•  • 

• 

+  .  .  .  + 

• 

• 

•  • 

X-j  (N)  ...  Xk(N) 

fK(°) 

X-j  (N-x) . . .  Xk(N-t) 

Vt) 

where 

0 ( t )  =  single  channel  output  trace 
X . ( t )  =  ith  input  channel 
K  =  number  of  input  channels 
f . (t)  =  filter  for  lag  t 

Equation  A. 1.1  may  be  expressed  as 


0(t)  = 

X1(l)f1(o)+...+X<(o)fK(o) 

X-|  (-t)^  (t)  +  .  .  .+XK(-i)fK(T) 

•  • 

+  .  .  .  + 

•  • 

•  • 

X1(N)f1(o)+...+XK(N)fK(o) 

•  • 

X] ( N-t ) f i (t)+. . . +XK( N-t ) f k(t ) 

X  ^ ( o ) f ^ (o)  +  . . .+X 


X-J  (Njf-J  (o)  +  . .  .+X 


Rewriting  equation  A. 
O(t)  =  X]  ( o ) f  1  ( o )  + . . .  +X-j 


( N )  f -j  ( o )  + . . .  +X-j 


This  may  be  expressed 
0(t)  =  X-j  ( o )  . . .  X1  (-t) 


X] (N)  ...  X] (N-t) 


1  38 


(o)fK(o)  +  . . .+X1  (-t)^ (t)  +  . . ,+XK(-x)fK(T) 


(N)fK(o)  +  . . .+X1  ( N-t ) f1 (t)  +  . . .+XK(N-T)fK(T) 

A  .1  .2 


.  2  ,  we  obtain 


( -T ) f 1 (t)+. . .+XK(o)fK(o)+. . . +XK( -t) f k(t) 


(N-x)f^ (t)+. . . +XK( N) f K(o )+. . .+XK(N-i)fK(T) 

A  .1  .3 


a  s 


f  -j  ( o ) 

• 

+ . .  .  + 

XK(o)  ...  xk(-t) 

•  • 

•  • 

Vo) 

• 

f-,  (t) 

•  • 

Xk(N)  ...  Xk(N-t) 

fK(x) 

A  .1  .4 


APPENDIX  2 

THREE  CHANNEL  MAXIMUM  LIKELIHOOD  FILTER 

The  following  discussion  of  the  three  channel 
maximum  likelihood  filter  is  based  largely  on  two  references 
(Mercado,  1973;  1978).  The  author  has  simply  tried  to 
clarify  some  of  the  steps  in  the  derivation. 

If  we  let  the  three  input  channels  be  denoted  by 
X  -|  ( t )  ,  X  ^  ( t )  and  X  ^  ( t )  ,  then  the  single  channel  output 
0  ( t )  is  the  convolution  of  the  input  traces  with  the 
filters  f-|(x),  f  £  (  t  )  and  f  3  ( x  )  .  Mathematically,  we  have 

o ( t )  =  X1 (t-x)*fi (x)+X2(t-T)*f2(T)+X3( t-l)*f3(T) 

A. 2.1 

Since  the  filter  is  to  pass  coherent  signal 
arriving  with  infinite  apparent  velocity,  the  filters  must 
satisfy 


f-|  (x)+f2(x)+f3(x)  =  5(x) 


A. 2. 2 


or 


f -j  (  x  )  =  6  (  x  )  -  f  2  (x  )  -  f  3  (  x  ) 


A  .2  .  3 


Substituting  equation  A. 2. 3.  into  A. 2.1 


‘ 


0(t)  =  X1(t-T)*[«(x)-f2(T)-f3(T)]+X2(t-T)*f2(T)+X3(t-T)*f3(T) 

=  X,  ( t-T ) *6 ( T )  —  X ,  (t-T)*f2(l)-X1  (t-T)*f3(T)+X2(t-T)*f2(t) 
+X3(t-x)*f3(x) 

=  X1  (t)+[X2(t-x)-X1 (t-x)]*f2(x)+[X3(t-x)-X1  (t-x)]*f3(x) 

A. 2. 4 

Now,  let 

Ai  (t-T)  =  Xi  (t-T)-X]  (t-x)  i  =  2,3  A. 2. 5 

Substituting  A. 2. 5  into  A. 2. 4, 

0(t)  =  X-j  (t)+A2(t-x)*f2(x)+A3(t-x)*f3(x)  A. 2. 6 

We  would  now  like  to  obtain  the  least-squares 
solution  for  f 2 ( t )  and  f  3  ( x  )  that  will  minimize  the  output 
power.  The  output  power  is  given  by 

P  =  0(t)T*0(t) 

=  [X] (t)+A2(t-x)*f2(x)+A3(t-x)*f3(x)]T*[X1  (t)  + 
A2(t-x)*f2(x)+A3(t-x)*f3(x)] 


s 
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=  [X1  (t)T  +  f2T(T)*A2(t-T)T  +  f3(T)T*A3(t-T)T]*[X1  (t)  + 
A2(t-T)*f2(T)+A3(t-T)*f3(T)] 

=  X1  (t)T*Xl  ( t )  +X.|  (t)T*A2(t-T)*f2(T)+X1T(T)*A3(t-T)*f3(x) 

+f2T(T)*A2(t-T)T*Xi(t)+f2T(T)*A2(t-x ) T*A2 ( t-x ) *f 2 (t ) 

+f2T(T)*A2(t-T)T*A3(t-T)*f3  (T)+f3T(T)*A3(t-T)T*X1 (t) 

+  f3T(x)*A3(t-T)T*A2(t-T)*f2(T  )+f3T(x)*A3(t-x)T*A3(t-x)*f3(x) 

A. 2. 7 

For  vectors  A  and  B,  and  matric  C;  the  product 
A*C*B  is  a  scalar,  and  thus,  A*C*B  =  b'*"*c"'"*a"'"  .  Using 
this,  equation  A. 2. 7  simplifies  to 

P  =  X]  (t)T*X1  (t)+2Xl  (t)T*A2(t-x)*f2(x)+2X1T(t)*A3(t-x)*f3(x) 

+f  2"*"  ( x )  *A2  ( t-x )  T*A2  ( t-x )  *f  2  (x )  +2f  2^(x )  *A2  ( t-x )  T*A3  ( t-x )  *f  3  (x ) 

+  f3"'"(x)*A3(t-x)"'"*A3(t-x)*f3(x)  A. 2. 8 

Using  arguments  analogous  to  those  found  in 
appendix  1,  we  have 

A2(t-x)*f2(x)+A3(t-x)*f3(x)  =  Af 


A. 2. 9 


* 


where  A  is  the  two  channel  convolution  matrix  of  A2 ( t ) 
and  A^ ( t ) 


and, 


A  = 


a2(o) 


a2(n) 


A3(°) 


A3(N) 


f 


f  2  (0  ) 
f3(0) 


f2  (t  ) 
f3(x) 


.  .  .  A2(-t)  A3(-t) 


A. 2. 10 


. .  .  A2 ( N-t  )  A^ ( N-t  ) 


Therefore,  equation  A. 2. 8  may  be  compactly 
written  as 


P  =  X1  (t)T*X1  (t)+2X1  (t)T*A*f  +  fT*AT*A*f  A. 2. 12 

The  minimization  of  the  power  and  the  least-squares 
solution  for  f  is  accomplished  by  differentiating  P  with 
respect  to  f ,  and  setting  the  result  equal  to  zero. 

3  (  X , (t)T*X,  (t)  2  8 ( X,  ( t ) T*A*f  9 ( f T* AT*A*  f ) 

ii-  =  - 1 - ! -  +  - ! -  +  - 

9 f  9  f  9  f  3f 

=  0  +  2X] (t)T*A+2fT*AT*A 


A. 2.1  3 


Therefore  , 


2  (  X  -j  ( t )  T*A  )  T  +  2  AT*A*f  =  0  A. 2. 14 

This  leads  to 

AT*A*f  =  -(X1  ( t )  T*A ) T  A. 2. 15 


that  is, 


f  =  -  (AT*A)_1  (X]  ( t ) T*A ) T  A. 2. 16 

At  this  point,  we  may  note  that  equation  A. 2. 16  is 
the  filter  equation  for  the  Wiener  filters.  The  main 
advantage  of  the  maximum  likelihood  filter  over  the  Wiener 
filter  is  that  for  N  input  channels,  N  filters  must  be 
calculated  for  Wiener  filtering  whereas  only  N-l  filters 
have  to  be  calculated  for  maximum  likelihood  filtering. 
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A3  .1 
A3. 2 
A3  .3 


APPENDIX  3 

LISTINGS  OF  COMPUTER  PROGRAMS 

The  velocity  filter 
The  n-th  root  stack 

The  three  channel  maximum  likelihood  filter 


A3.1  The  Velocity  Filter 

The  velocity  filter  program  used  in  this  thesis  was 
initially  written  by  Clowes  (1969)  as  part  of  his  Ph.D. 
thesis,  and  later  revised  by  the  author  for  analysis  of 
the  data  collected  in  southwestern  Manitoba.  This 
program  velocity  filters  seismic  data  according  to  the 
algorithm  located  in  section  4.1. 

For  use  in  the  program,  the  input  channels  can  have 
a  maximum  of  1500  points,  and  they  must  be  stored  in  trace 
mode  on  magnetic  tape.  This  input  data  is  read  off  of 
the  tape  through  the  use  of  the  system  routine  READ.  To 
plot  out  the  input  and  output  traces,  a  call  must  be  made 
to  the  subroutine  Plotvf,  which  makes  use  of  the  Calcomp 
plotter  subroutines. 


C  THE  VELOCITY  FILTER  PROGRAM 

REAL  *4  SUM  (3000)  ,Y(3000) ,Z  (3000) ,  TS  (64) ,RB  (1500) 
COMMON/S TOR A G/XDATA (15,3000) 

C  XDATA. . . THE  MATRIX  OF  INPUT  CHANNELS 

COMMON/STORE/XFILT (12,2500) ,SEPT, SC ALE, COPS , RECNO (5) 
1  , CENTER,  NTAU, NS, TSTART 
C  XFILT. . . THE  MATRIX  OF  OUTPUT  CHANNELS 

INTEGER  PTSTAR, ROWST, CENTER , START, CENT (5)  ,TAU(5) 

INTEGER*2  LEN 

XPOS=2 . 0 

YPOS=29. 0 

NTR  PT-0 

100  FORMAT  (1615) 

101  FORMAT (8F10. 0) 

103  FORM  AT (5 A4) 


READ (5, 100)  ISTOP 

C  ISTOP. ..IF  999, THE  TEST  PULSE  IS  USED. 

C  IF  100, THE  INPUT  TAPE  IS  USED. 

READ (5, 100)  M, NTR, NOUT, JPLOT, PTSTAR, NPEND, NS, NL 
C  JPLOT... IF  0 , NO  PLOT  IS  DONE 

C  IF  1, INPUT  AND  OUTPUT  TRACES  ARE  PLOTTED. 

C  NTR . NUMBER  OF  INPUT  TRACES 

C  M . NUMBER  OF  TRACES  TO  BE  USED  IN  THE  FAN  PASS 

C  OPERATION 

C  (M  MUST  BE  EVEN) 


C  NOUT. ...  NUMBER  OF  OUTPUT  TRACES 

READ (5,  101)  SEPT, COPS, SCALE, TSTART 
C  TSTART. .. START  TIME  FOR  AXIS  IN  PLOT 

C  SCALE. ...  SCALING  FACTOR  TO  BE  USED  IN  THE  PLOT 

C  SEPT . SEPARATION  BETWEEN  TRACES  IN  THE  PLOT 

C  COPS . NUMBER  OF  POINTS  PER  SECOND. 

READ  (5,  10  0)  (CENT  (I)  ,1  =  1,  NS)  ,  (TAU  (I)  ,1=1, NS) 

C  CENT  (I) _ THE  MOVEOUT  IN  POINTS  PEE  TRACE  ABOUT 

C  WHICH  THE  PASSBAND  IS  CENTERED. 

C  TAU (I) . . . . THE  MOVEOUT  IN  POINTS  PER  TRACE  OF  THE 

C  CUTOFF  VELOCITY 

C  NS . NUMBER  OF  PASSBANDS  TO  EXECUTE 

IF  (JPLOT. EQ.1)  CALL  PLOTS 
IF  (ISTOP .EQ. 999)  GO  TO  130 
NPTS= (NPEND-PTSTAR)  +1 

C  NPTS. .. NUMBER  OF  POINTS  IN  THE  INPUT  TRACES 

READ(5, 103)  RECNO 
C  RECNO. .. TITLE  FOR  THE  PLOT 

C 

C  READ  IN  THE  DATA  WITH  PTSTAR=1 
C 

DO  8  1=1 , NTR 

CALL  READ  (RB,LEN, 0,LX, 8) 

DO  8  J= 1 , NPTS 


‘ 

■ 


8  XDATA  (I ,  J)  =RB  (J) 

IF  (ISTOP.EQ.  100)  GO  TO  150 
C 

C  DO  THE  TEST  FUNCTION 
C 

130  CALL  PULSE  (TS) 

NPTS=NPEND-PTSTAR+ 1 

T  ST ART=0 • 0 

READ (5, 103)  RECNO 

153  IF  (JPLOT. EQ. 1) CALL  PLOTVF (NTR, NPTS, 1 , NOUT ,PTST'AR , 0 , 
1 NTR  PT,  XPOS,YPOS) 

C  PLOT  THE  INPUT  TRACES 
DO.  65  KF  =  1  ,NS 
CENTER=CENT (KF) 

N  TA U=TA  U  (KF) 

L  =  M-1 

MC=L*IABS  (CENTER) 

N EN D=N PEND -  (L+1)  *NTAU/2 
START=PTSTAR+ (L- 1) *NTAU/2 
IF  (CENTER. LT.0)  START= START+MC 
IF  (CENTER. GT.  0)  N END=N END-MC 
C  THE  STARTING  AND  ENDING  TIMES  HAVE  BEEN  SHIFTED  SO 
C  THAT  ONLY  POINTS  BETWEEN  PTSTAR  AND  NPFND  WILL  BE 
C  USED  WHEN  THE  DATA  POINTS  ARE  TIME  SHIFTED. 

A=0. 65465 
B=0. 98612 
C=  0.  13091 
D=  0.  2026424 
ROWST=0 

20  ROWST=ROWST+ 1 
DO  10  1=1, NL 
10  SUM  (I)  =0.0 

DO  50  K=START,NEND 

N  SN=K-PTSTAR+ 1 

I ROW=RO W  ST 

DO  40  J=  1 , M 

JA= J- 1 

MU=2*JA-L 

MUDE= (MU-1) *NTAU/2 
M  UAD= (MU+  1 ) *NTAU/2 
JB=CENTER*JA 
ICOL1=K-MUDE+JB 
ICOL2=K  +  MU  AD+JB 

DIFF= (XDATA (IRO W  ,ICOL 1 ) - XDATA  (IROW,ICOL2) ) /MU 
SUM  (NSN)  =SUM  (NSN)  +DIFF 
40  I ROW=IRO W+ 1 
50  CONTINUE 

NST1  =  STA  RT-PTST AR+ 1 
NEN 1=NEND-PTSTAF+1 

C  SUM  HAS  BEEN  FILLED  WITH  THE  WEIGHTED  SUM  OF  THE 


‘ 


n  o 


TIME  SHIFTED  TRACES 

DO  55  1=1, NL 
Y (I) =0. 0 
Z (I) =0.0 

55  XFILT (ROUST, I) =0.0 
DO  60  K=NST1 , NEN 1 
N  A=K— 1 
NB=K-2 

Y  (K)  =SUM  (NA)  -A*SUM  (NB)  +  B*Y  (NA)  -C*Y  (NB) 

NN=NEN1 -K+NST1 
NC=NN+ 1 
ND=NN+  2 

60  Z  (NN)  =  SUM  (NN )  -A*SUM  (NC)  +  B*Z  (NC)  -C*Z  (ND) 

DO  70  1=  NST 1 , NEN 1 
70  XFILT (RONST , I) =D*  (Z  (I)  -Y  (I)  ) 

C  THE  TRACES  HAVE  BEEN  FILTERED. 

IF  (ROEST.  LT.NOUT)  GO  TO  20 

IF  (JPLOT.EQ.1)  CALL  PLOTVF (ROWST, NP TS , 2 , NOUT , 1 , KF , 
1NTR  PT ,  XPOS ,  YPOS) 

C  PLOT  THE  OUTPUT  TRACES 
65  CONTINUE 

IF  (JPLOT.EQ.O)  GO  TO  120 
CALL  PLOT  (XPOS, YPOS- 28. 0,-3) 

CALL  PLOT  (0.  0,0.  0,999) 

120  STOP 
END 


SUBROUTINE  PLOTVF  (NTP , NPTS , I NDEX , NOUT, NP ST AR  , KF , 

; NTR  PT, XPOS, YPOS) 

C  OM MON/S TO  RAG/XDATA (1  5,3000) 

COMHON/STORE/XFILT (12,2500)  ,SEPT, SCALE , COPS , RECNO  (5) 
1 , CENTER, NTAU  ,NS ,T START 
INTEGER  CENTER 
NPEND=NPSTAR+NPTS-1 
T=2. 0/COPS 
XLEN=FLOAT  (NPTS)  *T 
A  BSINT=0 • 5 
MNTR=MOD  (NTR  PT,2) 

GO  TO  (20,30) , INDEX 
20  DO  27  J= 1 , NTP 

DO  22  JA=NPSTAR,NPEND 
22  XDATA(J,JA)=XDATA  (J,JA)  /SCALE 
27  CONTINUE 

CALL  PLOT  (XPOS, YPOS , -3) 

CALL  PLOT  (0.0, -1.0, -3) 

CALL  AXIS  (0. 0,0. 0, ' TIME  IN  SECONDS ',+ 1 5 , XLEN ,0.  0 , 


IT  START, ABSINT,1 0.0) 

YD=3.5+5. 5*S  EPT 

CALL  SYMBOL  (- 1 . 0 ,-YD ,0. 2 ,RECNO , 90 . 0 , 20) 

YDADD= YD+0 .50 

CALL  SYMBOL(-0. 5,-YDADD,0. 15  ,* DISTANCE  BETWEEN  TRACES 
1  IS  .146  KM.  1  ,90.0,35) 

GO  TO  39 

30  DO  10  J=  1 , NTP 

DO  5  JA=NPST AR, NPEND 
5  XFILT  (J,  JA)  =  XFI LT  ( J ,  J A)  /SCALE 
10  CONTINUE 

IF  (MNTR. NE. 0)  GO  TO  36 
33  CALL  PLOT  (- XLEN , 0. 0  , -3) 

36  IF  (KF.NE.1)  GO  TO  39 
CALL  PLOT  (0. 0,-2.  0,-3) 

CALL  AXIS  (0. 0,0. 0, ' TIME  IN  SECONDS ',+ 1 5 , XLEN ,0.  0 , 

1 TST  ART , ABSINT, 10.0) 

CALL  SYMBOL  (- 1 . 0 ,-0. 7  ,. 20 , ' VELOCITY ' , 90 . 0 , 8) 

CALL  SYMBOL  (- 0.  7  ,  - 0 .  7 , 0 . 2  0  ,  '  FI  LTERED '  ,9  0 . 0 , 8) 

CALL  SYMBOL  (-.  3,  - .  4,  0.  2  0,  1  DATA  »  ,  90. 0 , 4) 

39  IF  (KF.EQ.0)  GO  TO  44 

Y D=  2.  1 +  FLOAT  (NOUT-1)  *SEPT 
FNUM 1=  (  (CENTER-NTAU)  /COPS) 

IF  (FNUM1.NE.0.)  GO  TO  40 

Y  EL  1=  1 . 0 
GO  TO  41 

40  V EL 1=0. 293/FNUM 1 

41  FNUM2=  (  (CENTER  +  NTAU)  /COPS) 

IF  (FNUM  2. NE. 0. )  GO  TO  42 
V EL  2=1. 0 

GO  TO  43 

42  VEL2=0. 293/FNUM2 

43  CALL  SYMBOL(-1. 0, -YD,. 15, 'VI' ,90.0,2) 

CALL  SYMBOL  (-0.  9 ,-YD+ . 25 ,. 1 0 , ' APP '  ,9  0.0,3) 

CALL  WHERE  (XP,YP,FCTR) 

YP=YP+  0 • 2 

CALL  SYMBOL (-1. 0 ,YP,. 15, 1  IS  ',90.0,4) 

CALL  WHERE  (XP,YP,FCTR) 

Y  P=YP+  0.2 

CALL  NUMBER  (- 1 . 0 , YP , . 1 5 , VEL 1 , 90 . 0 , 1 ) 

CALL  WHERE  (XP,YP,FCTR) 

YP=YP+0 . 2 

CALL  SYMBOL  (-  1 . 0  ,  YP  ,  .  1 5  ,  '  K/S  •  ,  9  0 . 0 , 3) 

CALL  SYMBOL  (-0 .  5, -YD  , .  1  5,  '  V2  '  ,  90 . 0 , 2) 

CALL  SYMBOL  (-C  .  4, -YD  +  .  2  5,  .  1  0,  '  APP  ' ,  90. 0 , 3) 

CALL  WHERE  (XP,YP,FCTR) 

Y  P=YP+  0 . 2 

CALL  SYMBOL  (-0.  5,  YP  ,  .  1  5  ,  '  IS  ',90.0,4) 

CALL  WHERE  (XP,YP,FCTR) 

YP= YP+  0.2 


% 

- 


CALL  NUMBER  (-0 . 5 , YP , . 1 5 , VEL2  , 90 . 0 , 1 ) 
CALL  WHERE  (XP,YP,FCTR) 

YP=YP+0. 2 

CALL  SYMBOL  (-0. 5  ,  YP  , .  1  5  ,  •  K/S  '  ,  90. 0 , 3) 

44  DO  80  J=1,NTP 
M J  =  MOD  { J  ,  2) 

IF  (J.  NE.  1)  GO  TO  45 
CALL  PLOT  (0. 0,-2. 0,-3) 

GO  TO  46 

45  IF  (HJ.EQ.0)  GO  TO  60 
CALL  PLOT  (-XLEN  , -SEPT, -3) 

46  CALL  SYMBOL  (0. 0 , 0. 0 , . 15 ,04 , 0. 0 ,-4) 

GO  TO  (47,49) , INDEX 

47  CALL  PLOT  (0 . 0 ,  XD ATA  ( J ,  NPSTA R)  ,  3 ) 

DO  48  JA=N PS TAR, NPEND 

TSCALE  =  T*FLOAT ( J A- NPSTAR+  1 ) 

48  CALL  PLOT  (TSCALE , XDATA ( J , JA )  ,2) 

GO  TO  75 

49  CALL  PLOT  (0 . 0  ,  XFILT  ( J,  NPSTA  R)  ,  3) 

DO  50  JA=NPSTAR,NPEND 
TSCALE=T*FLO AT  (JA-NPSTAR+1) 

50  CALL  PLOT  (TSCALE  , XFILT  (J  ,JA)  ,2) 

GO  TO  7  5 

60  CALL  PLOT  (XLEN ,-SEPT ,-3 ) 

CALL  SYMBOL  (0. 0 , 0. 0 , . 1 5, 04 , 0. 0 , -4) 

GO  TO  (6  1,64)  , INDEX 

61  CALL  PLOT  (0. 0, XDATA (J,NPEND)  , 3) 

DO  62  JA=NPSTAR, NPEND 

T' SCALE=-T* FLOAT  (JA-NPSTAR+1) 

J  CC=NPEND— JA  +  NPSTAR- 1 
IF  (JCC. LT.NPSTAR)  JCC=NPSTAR 

62  CALL  PLOT  (T SCALE , XDATA (J , JCC) , 2) 

GO  TO  75 

64  CALL  PLOT  (0 . 0 ,  XFILT  (  J  ,  NPEND )  ,  3) 

DO  65  JA=NPSTAR, NPEND 

T  SCALE=- T*FLO AT (JA-NPSTAR+1) 

JCC=NPEND -JA+NPSTAR- 1 
IF  (JCC. LT.NPSTAR)  JCC=NPSTAR 

65  CALL  PLOT  (TSCALE, XFILT (J, JCC)  ,2) 

75  CALL  PLOT  (TSCALE ,0. 0 ,3) 

CALL  SYMBOL  (TSC ALE , 0 . 0 , . 1 5, 0 4 , 0 . 0 , - 4) 

80  CONTINUE 

GO  TO  (81,84) ,  INDEX 

81  DO  83  J= 1 , NTP 

DO  82  JA=NPSTAR, NPEND 

82  XDATA (J,JA) = XDATA (J, JA) * SCALE 

83  CONTINUE 
GO  TO  92 

84  DO  90  J  =  1 ,  NT P 

DO  85  JA=NPSTAR, NPEND 


' 


■ 


151 


85  XFILT  ( J ,  J  A) =XFILT (J, JA) * SCALE 
90  CONTINUE 

92  IF  (KF.NE.NS)  GO  TO  100 
IF  (MNTR.EQ.O)  GO  TO  95 
XPOS=XLE  N  +  5 . 0 
GO  TO  97 
95  XPOS=5 . 0 

97  YPOS=11. O  +  FLOAT  (3*NTP  +  8)  *SEPT 
100  NTS  PT=N TP 
RETURN 
END 


SUBROUTINE  PULSE  (TS) 

COMMON/S TOR AG/ XDATA (15,3000) 

REAL  TS (64) 

PI=3. 1415926 
DO  6  1=1,16 
JR=I-1 

TS  (I) =  -0  . 5*SIN (PI/1 6* JR) 

6  CONTINUE 

DO  7  1=17,32 
JP=I-17 

TS  (I) =SIN (PI/1 6* JP) 

7  CONTINUE 

DO  8  1=33,48 
JQ=I-33 

TS  (I) =-0  . 75*SIN (PI/16* JQ) 

8  CONTINUE 

DO  9  1=49,64 
JC=I-49 

TS  (I) =0. 25*SIN (PI/16* JC) 

9  CONTINUE 

C  DC  1  1=1,64 

C  JB=I-1 

C  TS  (I) =EXP (-PI/16 . *JR) *SIN (0.5*PI*JR) 

C  1  CONTINUE 

DO  2  J=1 ,  1  1 
DO  2  1=1 ,3000 
XDATA ( J, I) =0.0 
L  =  0 

DO  3  1=1,7 
11=1-1 

IB=  (200*1)  +1+  (11*64) 

L=L+  1 

DO  4  J= 1,11 
JA=  (L-4 )  *2 
JS=J-1 


2 


' 


JN=JS*JA 
IT=IB-JN 
NN=II+63 
DO  5  JT=IT  #NN 
JL=  (JT-IT)  +1 
5  XDATA (J, JT)=TS (JL) 
4  CONTINUE 
3  CONTINUE 
RETURN 
END 


A3. 2.  The  N-th  Root  Stack 

This  program  computes  the  eight  root  stack  for  a 
given  set  of  input  seismic  traces  according  to  the  theory 
found  in  section  4.2.  A  maximum  of  15  input  channels  is 
allowed,  and  each  trace  must  comprise  1500  data  points. 
These  input  traces  must  be  located  on  a  magnetic  tape 
with  a  block  length  of  6000  bytes,  and  are  read  with  the 
aid  of  the  system  routine  READ.  To  accommodate  different 
stepouts  in  the  stack,  the  first  and  last  50  points  of 
the  output  trace  are  set  equal  to  zero.  The  resulting 
output  traces  are  then  written  onto  magnetic  tape  via  the 


system  routine  WRITE. 


‘ 


non  non  non 


100 


THE  8TH -ROOT  STACKING  PROGRAM 
INTEGER  *2  L  EN , LEB 

REAL  *4  RB  (1  500)  ,Y  (1500)  ,YY  (1  500)  ,X(15 ,1500) 
RB...A  SEISMIC  CHANNEL  FROM  THE  INPOT  TAPE 
X....THE  MATRIX  OF  INPUT  CHANNELS 
YY...AN  OUTPUT  CHANNEL  TO  BE  WRITTEN  ONTO  TAPE 
FORMAT  (1  615) 

READ  (5,  10  0)  M, NTR, NS 

M.  .  .  .  0 NUMBER  OF  TRACES  TO  BE  USED  IN  THE  STACK. 

INPUT  TRACES 
POSITIVE  STEPOUTS 


OF 

OF 


READ 


NTR. . . NUMBER 
NS-1.  .  NUMBER 
LEB=60  00 
N  SS= (2*NS) -1 
ML= -NS 

NOTJT=  (NTR-M)  +1 
NOUT. .NUMBER  OF 
IN  THE  DATA. 


OUTPUT  TRACES 


DO  1  1=1 ,  NTR 

CALL  READ(RB,LEN,0,LX,8) 

DO  1  J=1  , 1  50  0 
1  X  (I,J)=RB(J) 

C  DO  NSS  EIGTH  ROOT  STACKS. 

DO  8  IR=1,NSS 
N  1=1 
N  2=M 
ML=ML+ 1 
DO  2  1=1, NOUT 
LQ=  50 

DO  3  J=1 ,  1500 

3  Y (J) =0. 0 
NN  =  0 

DO  4  JQ  =  N 1 ,N 2 
NN=NN+ 1 

DO  5  JC=5 1 ,1450 
LQ=  LQ+1 

IF  (X(JQ,LQ)  •  EQ.0. 0)  GO  TO  5 

Y  ( JC)  =  Y  (JC)  +  (ABS  (X  ( JQ  , LQ)  )  **0.  12  5)  *X  (JQ,LQ)  / 

;  ( AB S  (X  ( JQ  , LQ )  )  ) 

5  CONTINUE 

LQ=  (ML*NN)  +50 

4  CONTINUE 

DO  6  JC=1, 1500 

IF(Y(JC)  .EQ.  0.0)  YY  ( JC)  =0.0 

IF  (YY (JC) .EQ. 0.  0)  GO  TO  6 

YY  (JC)  =  (ABS  ( Y  (JC)  /FLOAT  (M)  )  **8)  *Y  ( JC)  /  (ABS  (Y  (JC)  )  ) 

6  CONTINUE 

CALL  WRITE  (YY, LEB, 0,LY,1) 

C  WRITE  THE  OUTPUT  ON  TAPE 

WRITE (6,100)  N1,N2,LQ,ML 


‘ 


- 
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N 1=N1+ 1 
N2=N2+  1 
2  CONTINUE 
8  CONTINUE 
STOP 
END 
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A3. 3.  The  Three  Channel  Maximum  Likelihood  Filter 

This  program  computes  a  three  channel  maximum 
likelihood  filter  according  to  the  theory  found  in 
appendix  2  .  A  maximum  of  eleven  input  channels  may  be 
used,  each  with  a  block  length  of  1000  bytes.  The  input 
traces  must  be  stored  in  trace  mode  on  magnetic  tape,  and 
are  read  with  the  system  routine  READ.  All  output  traces 
are  written  onto  tape  through  the  routine  WRITE. 

To  obtain  the  filter  coefficients,  equation  4.3.12 
must  be  solved,  which  is  of  the  form 

AX  =  B  A3. 3.1 


The  solution  to  this  set  of  linear  equations  is 
obtained  through  the  IMSL  library  routine  LEQTIF.  The 
subroutine  is  called  by 


CALL  LEQTIF(A,M,N,N,B,IDGT,WKAREA,IER) 

where 

A  =  input  autocorrelation  matrix  of  dimension  N  x  n. 

B  =  input  cross  correlation  matrix  of  dimension  N  x  1. 
On  output,  the  solution  X  replaces  B. 

WKAREA  =  work  area  of  dimension  N. 

IDGT  =  decimal  places  of  accuracy  in  the  elements 
of  A  and  B . 


IER  =  an  error  parameter 


Two  other  IMSL  library  routines  used  in  the  program  are 
VMULFP  and  VMULFM.  VMULEP  computes  the  matrix  product 
A*B^  through  the  call 

CALL  VMULFP  ( A  ,  B  ,  L  ,M  ,  N , L , N , C , L , I ER) 

where 

A  =  input  matrix  of  dimension  L  x  M 

B  =  input  matrix  of  dimension  N  x  M 

C  =  output  matrix  of  dimension  L  x  N 
IER  =  error  parameter 

VMULFM  computes  the  matrix  product  A^*B  by 
CALL  VMULFM  (A ,B , L ,M ,N , L,C,M  ,  IER) 

where 

A  =  input  matrix  of  dimension  L  x  M 

B  =  input  matrix  of  dimension  L  x  N 

C  =  output  matrix  of  dimension  M  x  N 
IER  =  error  parameter 

Other  subroutines  used  in  this  program  are  Trans 
(transpose  a  matrix).  Macron,  Multi  and  Plit  (plotting 
routine)  . 


' 


on  no  o  ooooooooonnoo  on 
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C  THE  3-CHANNEL  MAXIMUM  LIKELIHOOD  FILTER 

REAL  *4  DQ(40, 250)  ,DELT(2, 250)  ,RB(250)  ,B  (80,1)  , 

;  XX  (250)  ,  DEL  1  (2,  125)  ,  DEL2  (2,125)  ,  WKAREA  (8  0 )  ,  A  (80 , 80)  , 
;  XXX  (2,2  50)  ,TS(64)  ,R  (15)  ,X1  (12  5)  ,C1  (2,2)  ,  OTT  (250)  , 

;G1  (2,1)  ,G2  (2,1)  ,FF(40,1)  ,C2(2,2)  ,C22  (2,2)  ,CT  (250,1) 
INTEGER  *2  L EN ,LEB 
COM MON/STORAG/X  (11,250) 

COMMON/STORE/OT (9,250)  , N , NOUT , RECORD (5)  ,  RY (250) 

X.,.  ..THE  MATRIX  OF  INPUT  CHANNELS 
OT.. .THE  MATRIX  OF  OUTPUT  CHANNELS 

100  FORMAT  ( 1  61  5) 

101  FORMAT (4F10.  0) 

102  FORMAT  (6F15.  5) 

103  FORMAT  (' 0* , 10X, • FILTER  COEFFS.  ARE:  '/) 

111  FORMAT (5A4) 

211  FORMAT  ('  0*  ,1  OX,  5A4/) 

212  FORMAT('0',10X, 'SCALE  CHOSEN  ',F8.1/) 

R  E AD (5 , 10  0)  INDEX, JPLOT 
INDEX. ..IF  100, USE  THE  TEST  PULSE 
JPLOT.. .IF  100, A  PLOT  IS  MADE 

READ  (5, 100)  N,NC,NS,NF,NS1 ,NF1,NLAG,ISE, NSQ, NFQ 
N.... NUMBER  OF  INPUT  CHANNELS 
NC...3  FOR  A  3-CHANNEL  FILTER 
NS. ..DIGITAL  STARTING  POINT  FOR  EACH  TRACE 
NF... DIGITAL  END  POINT  FOR  EACH  TRACE 

NS1.. DIGITAL  STARTING  POINT  FOR  DESIGNING  THE  FILTER 
NF1.. DIGITAL  END  POINT  FOR  DESIGNING  THE  FILTER 
R  ZAD (5 , 101)  TLN , C  SEP , SCL , HT 
TLN... LENGTH  OF  THE  PLOT 

CSEP.. SEPARATION  BETWEEN  CHANNELS  IN  THE  PLOT 
SCL. ..SCALE  FACTOR  USED  IN  THE  PLOT 
HT.... HEIGHT  FACTOR  USED  IN  THE  PLOT 
READ  (5,  111)  RECORD 
RECORD. ..TITLE  OF  THE  PLOT 
LEB= 1000 
NOFS=NLAG+ 1 
N  B=NC- 1 
NOFS2=NB*NOFS 

CALCULATE  THE  NUMBER  OF  POINTS  PER  TRACE 
N  TS=  (NF-NS)  +  1 
NINT=N  TS- 1 
NTT= (NF1-NS1)  +1 

CALCULATE  THE  NUMBER  OF  POINTS  USED  TO  THE  DETERMINE  THE 
FILTERS 

SUBLN=TLN/  (FLOAT (NINT) ) 

NOUT=N-N  B 

NOFS... NUMBER  OF  OPERATOR  POINTS 
NOUT.. . NUMBER  OF  OUTPUT  TRACES 
X  POS=2. 0 
YPOS=29. 0 


* 


. 


n  n 
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BS2=- 1 . 0 
DO  91  J J= 1 , N OUT 
DO  91  J  =  1 , NT S 
91  OT(JJ,J)=0.0 
C  CALCULATE  RY  FOR  THE  PLOT 
DO  20  1=1, NTS 
20  RY  (I)  =FLOAT  (1-1)  *SUBLN 

IF (JPLOT. EQ. 100)  CALL  PLOTS 
IF  (JPLOT.  EQ.  100)  CALL  FACTOR(.8) 

IF  (INDEX. EQ. 100)  GO  TO  99 
C  READ  IN  THE  DATA 
DO  2  IQ=  1 , N 

CALL  READ  (RB,LEN,  0,LX,  8) 

HM=0 

DO  2  J Q=  NS , N  F 
MM=MM+1 

2  X  (IQ, EM)  =RB  ( JQ) 

GO  TO  30 

99  CALL  PULSE  (TS,R,  ISE) 

30  IF  (JPLOT. EQ.  100)  CALL  PLIT (TLN r CSEP, SCL , HT , NTS , 

;  1  ,XPOS,  YPOS) 

PLOT  THE  INPUT  TRACES 

DO  1  1=1  ,NOUT 
DO  10  J=1,NLAG 
LL= J+ 1 

DO  10  I J=  1 , N LAG 
DQ  (LL ,  I J )  =0.  0 
10  CONTINUE 
IR  1=  1 
11=0 

NQQ=NOFS 
DO  3  J= 1 , NTS 
11=1+1 
111=1+2 
XX  (J)  =X  (I,  J) 

DELT  (1,  J)=X  (II,  J)  -X  (I,  J) 

DELT (2, J) =X (III ,J)-X  (I,J) 

3  CONTINUE 
N  Q=0 

DO  7  IB=NS1, NF1 
NQ=NQ+ 1 

7  XI  (NQ)  =XX  (IB) 

C  CALCULATE  R0  AND  GO  AND  PLACE  IN  THE  AUTOCORRELATION 
C  AND  CROSS  CORRELATION  MATRICES 
C 

DO  4  J=  1  , N B 
MN=0 

DO  4  JJ=  NS  1, NF1 
MN=MN+ 1 


- 


4  DELI  (J,MN)  =DELT  (J,  JJ) 

CALL  VMULFP (DELI, DELI ,NB, NTT, NB, NB , NB , C 1 , NB  ,  IEE) 

CALL  MACR0N(A,C1,NB, N0FS2 ,NQQ,IR1,1, NTT) 

CALL  VMULFP(DEL1,X1,NB,NTT, 1 , NB , 1,G1,NB,IER) 

CALL  MULTI  (B , G 1 , IT , NB , 1 ,NTT,BS2) 

NMN  =  1 

C  CALCULAIE  THE  OTHER  AUTOCORRELATION  AND  CROSS  CORRELATION 
C  TERMS  AND  PLACE  THEM  IN  THEIR  RESPECTIVE  MATRICES 
DO  5  J= 1 , NLA G 
NQQ=NQQ-  1 
I E  1=IR 1  +  N  B 
NF2=NF 1-NMN 
N  S2=NS 1 - NM  N 
DO  6  JQ= 1 , NB 
MN  =  0 

DO  6  J J=NS2, NF2 
MN=MN+  1 

6  DEL2  (JQ,MN)  =DELT  (JQ,JJ) 

CALL  VMULFP (DELI , DEL 2 , NB , NTT, NB, N B, NB, C2 , NB , IER) 

CALL  TRANS  (C2,NB,C22) 

CALL  MACRON ( A , C2 , NB , NOFS2 , NQQ , IR 1 , 1 , NTT) 

CALL  MACRON (A, C22, NB, NOFS2 , NQQ, 1,IR1,NTT) 

CALL  VMULFP (DEL2, X1,NB, NTT, 1,NB,1 ,G2,NB,IER) 

CALL  MULTI  (B,G2, IT, NB, 1 , NTT, BS2) 

NMN=NMN+  1 

5  CONTINUE 

C  CALCULATE  THE  FILTERS 

CALL  LEQT1F (A, 1 , NOFS2 , NOFS2 , B, 5, WKAR EA , I ER) 

WRITE  (6, 103) 

WRITE (6, 102)  (B (IB , 1) ,IB=1,NOFS2) 

W  RITE (6 , 100)  NB,NS,NF,NOFS, NTT 
C  CALCULATE  THE  OUTPUT 
DC  11  J= 1 , NB 
ML=-1 

DO  26  J J=1 ,NOFS 

ML=ML+ 1 

MS=0 

NTS  1=NTS-ML 
DO  26  IQ= 1 , NTS 1 
M  S=MS+ 1 
MT=MS+ML 

26  DQ(JJ,MT)=DELT(J,IQ) 

DO  36  JQ=1,NOFS 
JB=  (  (2*JQ)  +  (J-1)  )-1 
36  FF  (JQ,  1)  =B  (JB,  1) 

CALL  VMULFM (DQ ,F F , NOF S , NTS , 1 , NOFS , NOFS , CT , NTS , IER) 

DO  47  J Q= 1 , N TS 
47  XXX  (J  ,  JQ)  =CT  (JQ,  1) 

11  CONTINUE 

DO  50  IH=NSQ,NFQ 
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50  OT  (I,IH) =  XX  (IH) +  XXX  (1 ,IH) +  XXX (2 ,IH) 

1  CONTINUE 

IF  (INDEX. EQ.  100)  GO  TO  41 
DO  92  J T= 1 , NOUT 
DO  93  JQ=1 ,NTS 
93  OIT  (JQ)  =OT  (JT,  JQ) 

92  CALL  WRITE  (OTT  ,LEB  ,  0  ,  LY  ,  1) 

41  IF  (JPLOT.EQ.  100)  CALL  PLIT (TLN,CSEP, SCL, HT,NTS, 
;  2  ,  XPOS  ,  YPOS) 

C  PLOT  THE  OUTPUT 

WRITE (6,211)  RECORD 
WRITE  (6,212)  SCL 
IF  (JPLOT.NE.  100)  GO  TO  42 
CALL  PLOT  (0. 0,0.  0,999) 

42  STOP 
END 


SUBROUTINE  PULSE  (TS , R , ISE) 
COMMON/STORAG/X (1 1 ,250) 
REAL  TS  (64)  ,  R  (15) 

PI=3.  1415926 
DO  6  1=1,16 
JR=I-1 

T  S (I)  =- 0. 5*SIN (PI/16* JR) 

6  CONTINUE 

DO  7  1=17,32 
J  P=I-1 7 

TS(I)=SIN  (PI/1  6* JP) 

7  CONTINUE 

DO  8  1=33,48 
J  Q=I-33 

TS(I)=-0.75*SIN  (PI/1  6* JQ) 

8  CONTINUE 

DO  9  1=49,64 
JC=I-49 

TS  (I)  =0.  25*SIN  (PI/1 6*  JC) 

9  CONTINUE 
R (1) =0. 0 
R (2) =0. 0 

R  (3)  =0.  25 
R (4) =0. 5 
R  (5)  =0.25 
R (6) =0. 0 
R  (7)  =-0.  25 
R (8) =-0. 5 
R  (9)  =-0.75 
R  (10)  =-1.0 


, 
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R  (1 1)=-0.75 
R  (12)  =-0,.  5 
R  (13)  =-0.25 
R  ( 1  4)  =0.  0 
R  (15) =0. 0 
DO  1  1=1,3 
DO  1  J= 1 , 2  50 
1  X(I,J)=0.0 
IL=0 

DO  3  10=1,3 
DO  4  1=1,3 
DO  5  J= 1 , 64 
IJ=  (1-1)  *64+34+IL+J 
5  X (IQ,IJ) =TS (J) 

4  CONTINUE 
IL=IL+I SE 
3  CONTINUE 
DO  10  1=1,3 
DO  11  J= 1 , 250 

IE(J.GT.  66.  AND.  J.LT.  8  2)  X  (I  ,  J)  =X  (I , J)  +  R  ( J-6  6) 

I  F  ( J. GT. 130. AND. J.LT. 146)  X  (I,  J)  =X  (I ,  J) +R  (J-130) 
11  CONTINUE 
10  CONTINUE 
RETURN 
END 


SUBROUTINE  PLIT (TLN , CSEP , SCL , HT , NTS, 
;INDIX,XPGS,YPOS) 

COMMON/STORAG/X  (11,250) 

CCMMON/STORE/OT (9,250) ,N , NOUT , RECORD (5) ,RY (250) 

T=TLN/250. 0 

XLEN=FLOAT  (NTS)  *T 

ABSINT=1.0/TLN 

TSTART=6 • 0+ (100. 0/250.0) 

GO  TO  (20,30) ,INDIX 
20  DO  27  J=  1 ,  N 
DO  22  J A=  1 , N TS 
22  X  (J,  JA)  =X  (J,  JA)  *HT/SCL 
27  CONTINUE 

CALL  PLOT  (XPOS , YPOS ,- 3) 

CALL  PLOT  (0. 0,-1. 0,-3) 

CALL  SYMBOL (0.0 ,1.5, 0.2, RE CORD, 0.0,20) 

CALL  SYMBOL(0.0, 1.0,0.15, ' DISTANCE  BETWEEN  TRACES 
;  IS  0.146  KM.  *,0.0,35) 

CALL  AXIS(0. 0,0.0, ‘TIME  IN  SECONDS' ,+1 5, XLEN ,0. 0 , 
; T  START, AESINT) 

GO  TO  45 


- 
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30  DO  10  J=1,NOUT 
DO  5  JA=1,NTS 

5  OT  (J,  JA)  =  (CT  (J,  JA)  *HT/SCL)  *5.0 
10  CONTINUE 

45  GO  TO  (40,50) ,INDIX 

40  CALL  PLOT  (0. 0,-2. 0,-3) 

DO  41  J= 1 , N 

DO  42  JA= 1 , NTS 

42  CALL  PLOT  (RY  (JA)  ,X  (J, JA)  ,2) 

CALI  PLOT  (0. 0,-CSEP,-3) 

41  CONTINUE 
GO  TO  46 

50  CALL  PLOT  (0.  0,-2.  0,-3) 

CALL  AXIS  (0. 0,0. 0, ' TIME  IN  SECONDS + 1 5, XLEN ,0 . 0 , 
; T  START, A ESIN T) 

CALL  PLOT  (0. 0,-2. 0,-3) 

DC  43  J= 1 , NO  UT 
DO  44  JA  = 1 ,NTS 

44  CALL  PLOT  (EY  (JA)  ,OT  (J,JA)  ,2) 

CALL  PLOT  (0.  0,-CSEP,-3) 

43  CONTINUE 

46  RETURN 
END 


SUBROUTINE  MULTI  (V,RR ,IS , NRA , NCA , NTT,BSN ) 
DIMENSION  V  (256)  ,RR  (NRA, NCA) 

DO  1  L= 1 , NCA 
DO  2  J=  1 , NRA 
IS=IS+1 

2  V  (IS) = (RR  (J, L) /FLOAT (NTT) ) *BSN 

1  CONTINUE 
RETURN 
END 


SUBROUTINE  MACRON  (A ,R 1 , NRA , NOFS2 , NST , NBL 1 , NEL2, NTT) 
DIMENSION  A  ( NOFS 2 , NO FS2)  ,R1  (NRA, NRA) 

I Z=NBL 1 -NRA 
I S=NBL  2-  1 
DO  1  1=1, NST 
IZ=IZ+NRA 
DO  2  J= 1 , NRA 
IS=IS+ 1 
IP=IZ-1 
DO  2  L= 1 , NRA 
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I F=IP+ 1 

2  A  (IS, IP)  =R1  (J,L)  /FLOAT  (NTT) 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  TRANS  (C2 , NB ,C22) 
DIMENSION  C2  (2,  2)  ,C22  (2, 2) 
DO  1  1=1, NB 
DO  1  J=1,NB 
1  C22  (J,I)=C2(I,J) 

RETURN 

END 


‘ 


